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I.  INTRODUCTION 


TKis  report  presents  a  three-dimensional,  vector,  boundary  element, 
or  surface  integral  moment-method  formulation  for  scattering  inside 
waveguides  and  represents  a  fairly  logical  extension  of  two-dimensional 
work  done  by  S.  Kagami  and  I*  Fukai,1  and  also  work  by  E.  Tonye  and  H. 
Baud rand. 2  The  method  uses  mode  expansions  in  the  regular  waveguide  of 
each  port,  and  allowance  for  multiple  ports  incurs  no  additional  theo¬ 
retical  cost,  as  this  is  a  multiport  formulation*  Therefore,  at  least 
in  principle,  multimode  multiport  scattering  matrix  information  can  be 
obtained*  The  primary  restrictions  are  that  the  discontinuity  boundary 
be  perfectly  conducting  and  that  the  material  filling  the  discontinuity 
volume  be  homogeneous  and  Isotropic.  These  restrictions  may,  to  some 
extent,  be  lifted.1 

Other  numerical  schemes  that  produce  multimode  scattering  data  are 
mode  matching1*  *5  at  abrupt  transitions  from  one  waveguide  to  another, 
and  other  methods6  based  upon  Solymar's7  equations  for  tapering  transi¬ 
tions  between  waveguides.  These  are  excellent  tools  that  have  recently 
been  developed  to  the  point  of  having  practical  application.  Both 
methods  apply  only  to  two-port  transitions,  and  the  central  axis  of  the 
two  connected  guides  must  be  parallel.  A  boundary  element  formulation 
overcomes  these  restrictions,  albeit  at  an  even  more  intensive  computa¬ 
tional  cost. 

Another  approach  that  could  be  potentially  more  flexible  than  the 
one  described  in  this  report  would  be  based  on  finite  element  or  finite 
difference  schemes.  Initially,  we  thought  a  boundary  element  scheme 


dimensional  problems,  but  in  three  dimensions  it  is  not.  It  turns  out 
that  for  a  boundary  element  scheme,  one  needs  to  compute  and  store 
coupling  coefficients  from  every  surface  patch  to  every  other  surface 


m. 


patch.  This  means  the  storage  required  is  proportional  to  the  surface 
area  squared  or  the  enclosed  volume  to  the  four-thirds  power,  whereas  a 
corresponding  difference  approach  requires  storage  proportional  to  the 
volume.  Without  knowing  the  coefficients  of  proportionality  for  the  two 
different  cases,  nothing  more  can  be  said  by  way  of  comparison.  It  is 
important  to  realize  however  that  a  boundary  element  formulation  in 
three  dimensions  does  not  a  priori  require  less  storage  than  a  finite 
element  formulation. 

An  interesting  and  possibly  less  computationally  intensive 
approach,  based  on  variational  techniques  for  irregular  waveguide  tran¬ 
sitions,  is  presented  by  Bernstein,  et  al.8 

This  work  was  motivated  by  the  need  for  a  tool  capable  of  analyz¬ 
ing  the  rapid  transition  from  a  regular  waveguide  to  a  slow  wave  cir¬ 
cuit.  Although  the  theory  presented  is  not  sufficiently  well  developed 
to  do  this  job,  the  goal  was  simply  to  develop  a  framework  upon  which 


such  a  problem  could  properly  be  considered.  This  is  primarily  an 


exposition  on  the  theoretical  and  computational  aspects  of  the  boundary 


element  method  as  it  applies  to  waveguide  discontinuities.  Numerical 


results  are  given  for  a  few  simple  but  important  test  cases  to  verify 


that  the  work  given  here  does  have  substance  and  to  demonstrate,  rather 


specifically,  what  types  of  computational  expense  to  expect. 
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II.  DERIVATION  OF  AN  APPROPRIATE  SURFACE  INTEGRAL  EQUATION 


In  this  section.  Green’s  theorem  Is  applied  to  a  scalar  wave 
function,  thereby  obtaining  the  wave  function  interior  to  a  volume  in 
terms  of  its  boundary  values.  The  point  of  evaluation  is  then  " pushed” 
to  the  surface  of  the  volume  by  a  limiting  process,  yielding  a  surface 
integral  equation  for  the  wave  function.  It  is  then  noted  that  the 
electric  and  magnetic  field  Intensities  satisfy  vector  wave  equations, 
and  that  the  result  derived  for  a  scalar  wave  function  applies  to  the 
rectangular  components  of  these  equations.  Vector  integral  equations 
are  then  obtained  by  recombining  the  component  equations.  To  wrap  up, 
the  integrand  is  simplified  for  the  case  when  part  of  the  boundary  is 
perfectly  conducting. 

Green '8  theorem. 


^  (♦Vi|i  -  +V$)  •  nda  -  /  ($V  <|i  -  tpp  d  x 

S  V 


applies  on  a  volume  V  with  volume  element  d  x,  S  is  a  closed  surface 

a 

bounding  V,  with  area  element  da  and  unit  outward  normal  n  at  da. 
Take  <|»  to  be  a  solution  to  the  inhomogeneous  wave  equation. 


(V2  +  k2)<Kx)  -  -  F(x) 


for  all  points  x  in  the  volume  V.  Solving  for  V  4*  in  this  equation  and 
substituting  the  result  into  Eq.  1  gives 
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mt 


I  _ Si /•-  «  -N  A  ,  *,  v  '  «  v  .*  '  ■  »  .  V  V  . 


■*.  s. 


’  4  '  *  -  »  "  »  J  • 

*  .  *  .  *  * 

t  '«  ^4  •*■**.  *  .  *  ' 


£  {♦(x)V*(x)  -  ♦(x)V*(x)}  •  nda  -  -  /  {F(x)$(x)  +  »(x)[V2  +  k2]*(x)}  d3x 
S  V 

(3) 

Choice  of  the  function  $(x)  is  quite  arbitrary,  but  Eq.  3  indicates  that 
some  reduction  in  complexity  is  possible  by  choosing  +  such  that 

(V2  +  k2)*(x)  -  -  6(^  -  x’>  (4) 


where  $(x  -  x')  is  the  Dirac  delta  function.  In  effect,  take  $(x)  ■ 
+(x,  x')  to  be  a  "Green’s  function."  If  the  geometry  under  considera¬ 
tion  is  excessively  simple,  then  the  geometrically  correct  Green's 
function,  satisfying  either  Dlrelchlet  or  Neumann  boundary  conditions, 
is  known  explicitly.  For  this  case,  Eq.  3  reduces  to  the  well  known 
result. 


♦  <x’)  -  /  F(iH(x,  ;*)  d3x 
V 


(5) 


Usually,  the  physically  appropriate  Green's  function  is  beyond  knowing 
explicitly,  and  so  it  makes  sense  to  choose  $(x,  x')  such  that  it  has  no 
geometrical  bias.  Clearly,  the  "best"  choice  in  this  case  is  the  free 
space  Green's  function. 


♦(x,  x') 


etjk|^'-x 

4w|x'-x 


(6) 


In  making  this  choice,  it  has  been  assumed  that  k  is  constant  throughout 
V.  As  you  might  have  guessed,  the  ±  sign  in  Eq.  6  is  arbitrary.  But 


for  e  time  dependence,  the  -  sign  is  the  standard  choice  and  the  one 


used  here*  For  future  reference, 

-jkR  r 

V$(x,  x')  -  - — y~  +  ^  “  x)  (7) 

4xR  L  J 

where  R  «  |x'  -  x|  •  When  Eq.  6  is  used  in  Eq.  3  and  x'  is  taken  inside 
V,  the  result  is 

t|>(x')  -  /  F(x)<fr(x,  x')  d3x  +  ^  {$(x,  xf  )V4»(x)  -  i|»(x)V<Kx,  x')}  •  n(x)  do 
V  S 

(8) 

For  x*  outside  V,  the  left-hand  side  of  this  equation  is  zero.  Compar¬ 
ing  Eq.  8  to  the  result  when  the  geometrically  correct  Green's  function 
is  known  (Eq.  5),  we  see  that  the  price  paid  for  not  having  the  correct 
Green's  function  is  in  the  addition  of  a  surface  integral  to  the  inte¬ 
gral  equation. 

Next,  the  evaluation  point  x'  in  Eq.  8  is  specialized  to  a  point 
on  S  by  a  limiting  process  that  is  outlined  by  Brebbia.3  Figure  1  is  a 


A 

n 


Fig.  1.  Geometry  for  specializing  x'  to  a  smooth  surface  S. 
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blown-up  view  depicting  a  small  section  of  the  smooth  surface  S  dis¬ 


torted  outward  by  a  small  hemisphere  having  radius  e  and  surface  Sg . 

The  point  x'  is  taken  to  be  at  the  center  of  the  hemisphere  and  on  the 

original  surface  S.  That  part  of  the  original  surface  now  displaced 

by  S  will  be  referred  to  as  A  .  in  this  geometry,  x’  is  inside  the 
6  & 

distorted  volume,  so  that  Eq.  8  remains  valid.  Now  let  us  see  what 
change  the  terms  in  this  equation  go  through  in  the  limit  as  the  dis¬ 
torting  radius  e  goes  to  zero.  Consider  the  surface  integration  in  the 
far  right  term  of  Eq.  8, 


lim  /  <Kx)V$(x,  x')  •  nda  +  lim  /  <|*(x)V$(x>  x*)  •  nda 

6+o  S-A  e+o  S 

e  e 


2*  x/2  _  .  -jker  n 

■  f  +(x)V<Kx,  x')  *  +  li®  /  /  'Kx*  +  er)  - j- Hk  +  —  (-e) 

S  6+0  O  O  L.  J 


-  e-jke 
4xe 4 


•  e  sin  0  d0  d$ 


f  <Kx)V<|)(x,  x')  •  nda  -  ~  <Kx’) 
S 


The  1/2  i|»  term,  resulting  from  the  limiting  process  above,  is  the 
classic  result  obtained  when  x'  is  specialized  to  any  smooth  point  on 
S.  When  specializing  to  corners  or  edges  on  S,  a  more  general  term  of 
the  form  nty  is  obtained,  where  n  is  between  0  and  1.  The  explicit  value 
of  n  can  be  obtained  for  each  case  by  the  method  used  here  for  smooth 
surfaces.  For  edges,  we  have  n  ■  0/2x,  where  0  is  the  edge  angle  as 
measured  from  outside  the  enclosing  volume.  For  right  angle  corners, 
h  -  1/8  and  7/8  for  inside  and  outside  corners,  respectively.  Turning 
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our  attention  now  to  the  two  remaining  integrals  in  Eq.  8,  we  see  that 
the  limiting  process  produces  no  more  new  terms  because  the  order  of  the 
singularity  in  the  integrands  is  too  low  in  both  cases  •  With  these 
results ,  Eq.  8  can  be  specialized  to 


■i<Kx')  -  /  F(x)$(x,  x’)  d3x  +  f  {$(x,  x')V^(x)  -  i|>(x)V4>(x,  x')}  •  ndo 

V  S 

(9) 

where  x'  must  now  belong  only  to  the  smooth  subset  of  points  making  up 
S. 

It  was  mentioned  earlier  that  the  electric  and  magnetic  field 
intensities  satisfy  vector  wave  equations.  This  is  simple  to  show  and, 
in  doing  so,  an  important  assumption  is  made  obvious.  The  harmonic, 
e^°)t,  form  for  Maxwell's  equations  is 


7  x  E  *  -  jwpH 


7  x  H  -  J  +  jweE 


For  p  and  e  constant 


(10) 


(ID 


7  •  (eE)  -  p 


7  •  (pH)  -  0 


7  •  E  -  p/e 


7  •  H  -  0 


(12) 


(13) 


Take  the  curl  of  Eq.  11  and  assume  e  constant. 


VxVxh-V*J+  j»e7  x  E 


Substitute  for  7  x  g  from  Eq.  10, 


2- 

V  x  V  x  h  -  «  weH  -  7  x  j 


Use  the  vector  Identity, 


7  x  7  x  H  -  7(7  •  H)  -  7  H 


to  obtain 


7(7  •  H)  -  71 


2H  -  k2H  - 


7  x  j 


2  2 

where  k  :  u  pe<  Now  assume  constant  y  and  use  Eq.  13  to  get 


2-  2- 

7H  +  kH--7xj 


The  parallel  result  for  the  electric  field  is 


2-2-  -  1 
7  E  +  k  E  -  jwyJ  +  -  7p 


Equations  14  and  15  are  the  inhomogeneous,  vector,  wave  equations  for 
the  electric  and  magnetic  field  intensities,  and  they  are  valid  at  all 
points  where  y  and  e  are  constant,  i.e.,  homogeneous  and  isotropic 
regions. 
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All  that  remain*  to  show  now  i*  how  the  result  for  a  scalar  wave 


function,  Eq.  9,  can  be  applied  to  the  vector  wave  Eqs.  14  and  15 «  To 
do  so  depends  only  on  the  fact  that  in  a  rectangular  coordinate  system, 
the  Lapladan  operator  may  be  written 


V2H  -  e^Hj  ♦  «272h2  *  *3v2h3 


where  e^,  ej,  and  e^  form  a  right-handed  set  of  orthogonal  unit  vectors 
in  a  rectangular  system.  In  component  form,  we  may  then  write  Eq.  14  as 


V2H±  +  k2^  -  -  [9  x  jj 


where  i  -  1,  2,  and  3.  This  is  a  scalar  wave  equation  to  which  we  can 
Immediately  apply  Eq.  9  to  yield 


y  H^x')  -  /  [V  x  J(x)]1  9(x,  x’)  d3x 


+  f  4(x,  x’)  ~  H1(x)  -  (x )  4(x,  x')  da 

S 


Recombining  the  component  parts  of  this  equation,  it  is  not  difficult  to 
show  that,  in  general, 


\  HU’)  -  /  [V  x  j(x)]  ♦<£,  x’)  dJx 
V 


-  -  3 


+  f  {*(x,  x’)[n  •  Vj  H(x)  -  H(x)[n  •  V*(x,  x* ) J }  da  (16) 
S 


-  9  - 


W 


i»,V, 


wi'rwivwvPA'.i, 


Equation  16  la  a  coordinate  fraa  raault  despite  recourse  to  a  rectangu¬ 
lar  system  for  Its  derivation.  It  can  be  applied  In  any  coordinate 


system.  The  parallel  result  to  Eq.  16  for  the  electric  field  can  be 
written  down  by  Inspection. 

Up  until  this  point,  the  source  tens,  J  and  p,  in  Maxwell's 
equations  have  been  retained  throughout.  There  are  two  reasons  for 
this,  one  is  simply  that  the  added  generality  required  very  little 

additional  work  and,  secondly,  because  the  equations  aay  be  useful  for 
exaalnlng  cavity  or  circuit  interaction  with  an  electron  beam.  With 
regard  to  this  second  reason,  great  care  must  be  exercised  in  the  appli¬ 
cation  of  an  equation  such  as  Eq.  16  because,  in  high  space  charge 

regions,  the  effective  peraittivlty  nay  vary  with  poeitlon  and  direc¬ 
tion;  l.a.,  plasma  permittivities  are  usually  tensor  quantities.  So, 
except  in  low  space  charge  devices,  the  equations  of  this  section  aay 
not  be  applicable.  From  here  on,  it  will  be  assumed  that  there  are  no 
source  terms  in  V. 

The  only  thing  remaining  to  this  section  is  to  show  how  the  inte¬ 
grand  of  Eq.  16  simplifies  when  all  or  pert  of  the  bounding  surface,  S, 
is  perfectly  conducting.  At  the  surface  of  a  perfect  conductor,  we  know 
that  the  tangential  components  of  the  electric  field  and  the  normal 

component  of  the  magnetic  field  are  zero;  E_  *0,  H  >0.  So,  using  Eq. 

t  n 

11  with  j  -  0, 


V  x  H  «  jucE  n 
t  J  n 
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this,  it  follows  that  ths  normal  dsri votive  of  ths  tangential 


components  of  H  are  0;  5fc  •  0.  Lot  that  part  of  the  surface  that  is 
perfectly  conducting  be  represented  by  8e  and  represent  the  remaining 
pert  of  the  enclosing  surfaca  by  Sf,  so  that  symbolically  S  -  Sc  +  Sf. 
With  this  understanding,  Bq.  16  may  be  written  as 

t  5  -  *  (I:  '  5t  (U)l <•“  *  f  [♦  (I? 5)  - 5  (D:)  <17> 

c  l  J  f  l  J 

where  now  the  volume  integration  has  been  dropped  and  the  explicit 
evaluation  points  are  implied  by  reference  back  to  Eq.  16.  Note  that 
the  integral  over  Sc  contains  only  three  unknown  functions,  the  two 
transverse  components  of  H  and  the  normal  derivative  of  the  normal 
component  of  H,  as  opposed  to  what  appears  to  be  six  unknown  functions 
on  the  general  boundary  Sf . 

In  the  next  section,  some  approximations  will  be  made  regarding 
the  nature  of  the  unknown  functions  on  both  boundary  types,  and  the 
number  of  algebraic  equations  required  to  find  these  functions  will  be 

made  clear. 


Cv:* 
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III.  APPLICATION  OF  THE  METHOD  OP  MOMENTS 


The  Intent  of  this  section  is  to  demonstrate  a  method  for  extract* 
ing  useful  Information  from  the  vector  surface  integral  equation  of 
Section  11,  Eq.  17.  The  approach  outlined  here  is  just  one  very  simple 
application  of  the  general  technique  called  the  method  of  moments,  or 
HoM.9 

First,  the  field  boundaries  are  considered  to  be  cross  sections  of 
regular  uniform  waveguides  for  which  the  normal  node  functions  are 
known,  and  the  surface  connecting  the  uniform  waveguides,  the  Junction 
surface,  la  assumed  perfectly  conducting  in  keeping  with  Section  II. 
Pulse  basis  sets  are  used  to  expand  the  unknown  H  field  on  the  conduct¬ 
ing  boundaries,  and  normal  mode  expansions  are  uaed  on  the  field  bound- 
erles  in  the  regular  waveguides.  The  system  is  excited,  or  driven,  by  a 
pure  mode  incident  on  the  junction  from  any  one  of  the  porta.  An  over¬ 
determined  linear  system  of  equations  Is  obtained  by  delta  testing, 
point  sMtchlng,  on  both  the  conducting  and  field  boundaries.  The 
unknown  set  contains  the  scattered  mode  amplitudes  at  all  the  ports 
along  with  the  field  amplitudes  for  the  pulse  basis  expansions  on  the 
conducting  boundaries. 

Figure  2  illustrates  a  general  N  port  junction  and  Indicates  how 
the  discontinuity  volume  Is  defined  by  placement  of  planar  field  bound¬ 
aries  perpendicular  to  the  axis  of  each  regular  waveguide  entering  the 
junction.  These  field  boundaries  can  be  moved  far  enough  away  from  the 
junction  and  Into  the  waveguide  so  that  only  propagating  modes  need  to 
be  considered.  Alternatively,  they  can  be  moved  as  close  to  the 
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rig.  2*  Schematic  diagram  of  an  N  port  junction  Indicating  placement 
of  field  boundaries  defining  the  discontinuity  volume  V. 


junction  as  possible,  while  still  remaining  inside  the  waveguide,  by 
including  an  appropriately  large  set  of  evanescent  nodes  along  with  the 
propagating  nodes  on  the  field  plane.  This  minimises  the  discontinuity 


Explicitly,  formal  node  expanalons  can  be  used  to  express  the 
field  at  boundary  j. 


HJ  "  6iJhU.€J  +  j^n]*[J.nJh[j.nl 


where  6^  is  the  Kronecker  delta  function  and  is  used  to  indlcste  that 
the  incident  node  is  in  the  ith  waveguide.  Other  definitions  are: 


Indices  over  field  boundaries. 

Integer  indices  used  for  ordering  and  counting  modes. 
Each  index  corresponds  to  a  unique  node  at  each  field 
boundary.  Mode  {  at  field  plane  1  is  the  incident  or 
driving  node. 
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Nona  11  tad  mode  function  at  flald  boundary  1  and  for 
aoda  n.  Tha  back  arrow  ia  uaad  to  Indicate  that  thia 


hU,il 

aoda  la  propagating  In  tha  -  n  direction.  Thia  func¬ 
tion  haa  tha  unite  of  magnetic  flald  intanalty. 

♦ 

h,.  ,  Same  aa  above,  except  the  forward  arrow  indicates 

ll»l  1 

A 

propagation  in  the  +  n  direction.  This  function  is 
used  for  all  reflected  and  transmitted  nodes* 

a,.,  .  Amplitude  of  node  1  at  field  plane  i.  Carries  no 

(1,1  ] 

units,  but  the  power  carried  by  the  associated  mode  is 
proportional  to  aa*. 

A 

The  distinction  between  s»de  functions  propagating  along  +  or  -  n  direc¬ 
tions  in  the  regular  waveguides  is  necessary  due  to  a  subtle  difference 
between  them  as  a  result  of  assumed  propagation  direction.  This  dis¬ 
tinction  la  Important;  it  means  the  difference  between  success  and 
failure.  Appendix  A  clerifles  this  distinction  and  also  gives  details 
on  the  normalisation  of  the  mode  functions  both  above  and  below  cutoff. 

Equation  17  also  requires  the  normal  derivative  of  h  on  the  field 
boundaries.  To  this  end. 


a__ 

in 


■j4iJBlJ,*)hU.tl  "  J 


L  • 

U  ,i) 


(19) 


where  ft,.  .  is  the  propagation  constant  for  mode  1  at  field  plane  j  and 
IJ  ft*  J 

will  be  purely  imaginary  below  cutoff.  The  formula  for  coeputlng  the 
ft's  are  given  In  Appendix  A. 
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New  suppose  that  the  conducting  surface  can  be  broken  down  into  a 
set  of  staple  geometric  regions.  That  is,  subdivide  the  conducting 
surface  into  rectangles,  triangles,  sections  of  cylinders,  etc.,  such 
that  all  of  the  regions  combined  describe  or  approximate  the  conducting 
surface  of  interest.  The  point  of  this  is  to  use  shapes  that  can  be 
defined  on  a  computer  with  a  minimum  of  effort.  On  each  region,  define 
a  local  set  of  orthogonal  curvilinear  surface  coordinates  (p,s).  The 
local  coordinates  correspond  uniquely  to  a  point  In  the  absolute  coordi¬ 
nate  system  (x,y,z)  by  the  surface  mapping  M,  which  depends  on  the 
surface  shape,  orientation,  and  location.  Symbolically,  M(p,s)  ♦ 
(x,y,z)  for  all  p,a  on  the  defined  region;  i.e.,  M  takes  the  surface 
coordinates  (p,s)  into  the  absolute  system  (x,y,z).  The  surface  tangent 

A  *  A  A 

unit  vectors  should  be  such  that  p  *  s  -  n,  where  n  is  the  outward  unit 

A  A 

normal  to  the  enclosed  volume  V.  The  unit  vectors  p  and  s  are  later 
referred  to  as  the  unit  primary  and  unit  secondary  directions. 

The  transverse  components  of  H  on  the  conducting  surface  can  now 
be  expanded,  approximated,  by  a  pulse  basis  expansion  on  each  region. 


H  <a;p,s)  -  I  H  Pf  kn(p,s) 

C  [a.kl]  C [a,kl 1 


where  the  following  definitions  are  used: 


Integer  index  used  for  ordering  and  counting  the 
simple  geometric  regions  used  to  build  the  conducting 
surface. 


-  15  - 


VV’.»  V.V/ 


k,l  Patch  position  Indices  used  for  ordering  and  counting 

the  patch  positions  on  region  a* 

Local  coordinates  on  region  a. 

Sum  over  patch  positions  k,l  on  region  a. 

Constant,  transverse  magnetic  field  amplitude  at 
position  k(l  on  region  a.  Carries  units  of  aagnetic 
field  intensity. 

Hc(a;pts)  Approximation  to  transverse  H  field  at  local  coordi¬ 
nate  positions  p,s  on  region  a. 

P[a  1  if  P.s  is  on  patch  k,l  in  region  a.  0  otherwise. 

This  last  equation  defines  the  pulse  function.  Note  that  nothing  here 

requires  the  patch  shape  to  be  rectangular  or  flat;  indeed,  the  sub- 

doaaln  patches  ordered  by  Indices  k,l  on  region  a  need  not  be  flat  or 

rectangular.  The  pulse  function  is  simply  a  function  having  value  1 

when  the  local  coordinates  p,s  are  Inside  patch  k,l,  and  0  otherwise. 

The  normal  derivative  of  the  normal  component  of  H  can  also  be  expanded 

a 

on  a  pulse  basis.  Let  H'  =  t—  H  .  Then, 
r  n  an  n 


P»* 

i 

l*»kl] 

Ht 

[a.kl] 


H’(a;p,s)  -  l  H’  P,  V11(P»*> 

n  [a.klj  n[a,kl]  [a*kl] 


where  all  the  same  notation  applies. 


While  the  pulse  basis  set  may  be  the  most  flexible  and  convenient, 
it  is  certainly  not  the  most  accurate  basis  expansion  available.  The 


pulse  expansions  In  Eqs.  20  and  21  approximate  smooth  functions  by 
functions  with  abrupt  jumps  at  patch  boundaries.  Clearly,  some  physics 
has  been  lost  by  this  approximation.  Reasonably  simple  subdomain  basis 
sets  that  are  capable  of  C°  approximations,10  smooth  approximations,  are 
linear  interpolation  on  a  triangular  net  and  bilinear  interpolation  on  a 
rectangular  net.  On  a  triangular  patch,  the  linear  interpolating  poly¬ 
nomial  g  has  the  general  form,  g  -  a  +  bp  +  cs,  where  p  and  s  are  the 
surface  primary  and  secondary  coodinates  as  before.  The  three  coeffici¬ 
ents,  a,  b,  and  c,  are  uniquely  related  to  the  values  of  g  at  the  three 
vertices  of  the  triangle  by  the  index  equation  g^  -  a  +  bpA  +  csjj  1  - 
1,  2,  3.  This  is  trivial  to  invert  for  the  polynomial  coefficients  in 
terms  of  the  vertex  values  g^.  The  underlying  basis  function  for  this 
interpolant  is  sometimes  called  the  rooftop  function  and  can  be  found  by 
setting  two  of  the  three  values  for  g^  to  0  and  the  third  value  to  1, 
then  solving  for  a,  b,  and  c.  For  rectangular  patches,  a  bilinear 
interpolating  polynomial,  g  ■  a  +  bp  +  cs  +  dps,  is  used  and  again  the 
four  coefficients  are  uniquely  determined  by  the  four  vertex  values. 
The  underlying  basis  for  this  Interpolant  is  called  the  pagoda  function 
and  can  be  found  as  before  by  setting  three  of  the  four  vertex  values  to 
0,  while  setting  the  fourth  equal  to  1.  These  basis  functions  are 
almost  as  flexible  as  pulse  functions.  The  additional  organizational 
complexity  required  to  implement  a  scheme  using  more  appropriate  basis 
sets  may  prove  worthwhile  by  producing  more  accurate  solutions. 


At  this  point,  we  have  approximated  all  the  surface  fields  using 
expansion  functions.  This  limits  the  domain  of  the  function  space.  Now 
we  have  only  a  finite  set  of  unknown  coefficients  to  look  for,  not 


entire  functions*  On  the  field  boundaries ,  the  unknowns  are  the  mode 
amplitudes,  and  on  the  conducting  boundaries,  the  unknowns  are  the  pulse 
function  aaplltudes.  The  next  job  is  to  relate  these  unknowns  through 
the  physics  of  the  situation. 

Clearly,  the  surface  integral  equation  from  Section  II,  Eq.  17, 
provides  the  relationship  we  need  between  the  surface  fields.  Using  the 
field  expansions,  Eqs.  18-19  and  20-21,  in  the  right-hand  side  of  Eq.  17 
yields  a  formidable  looking  result 


i'V  j, 


i  H  -  l  l  '  h;  f  ♦  n  do 

a  [a,kl]  [a.kl]  [a.kl] 


-  H  j-  p  (7$  •  n)  da 

P[a,kl]  [a.kl] 


-  H  f  s  (V*  •  n)  da 

[a.kl]  [a.kl] 


+  fSi[J8[1’cl* ' v*  ’  do 
-  I  l  a[1  nl  t  '  jB[1  Til*  +  7*  *  M 

J  Ij.n]  lJ,T1J  s  u*nJ 

J  V  ^ 


[j  »n] 


Where  the  transverse  H  field  has  been  decomposed  into  its  primary  and 
secondary  surface  components. 


H  «  pH  +  sH 

t[a,kl]  P[a,kl]  ®[a,kl] 


Note  that  the  integrands  contain  only  explicitly  known  functions;  the 
free  space  Green's  function  is  given  by  Eq.  6  and  its  gradient  by  Eq.  7. 
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We  can  now  generate  a  systems  of  linear  equations  by  systemati¬ 
cally  enforcing  Eq.  22  for  a  discrete  number  of  points  in  the  range 
of  x'  on  S.  This  is  called  delta  testing  or  point  matching.  So,  test 
Eq.  22  at  the  center  of  every  patch  on  each  region,  points  [c,mn],  and 
use  expansion  Eq.  20  on  the  left-hand  side  to  obtain  three  scalar  equa¬ 
tions,  one  for  each  of  the  three  vector  components.  They  are: 

The  primary  equation. 


-[c,mn]  n  r  *  “1C, 

p[c,mn]  *  li,£]  “  j  j®[j  ,n]P[c,mn]  *  F[j  , 


-[c,mn] 
n] 


+  l  l 

a  [a,kl] 


.  g[c,mn] 


’  H”[a,kl)Ptc’mnl  ’  ^a’klj 


.  H  *  .  -fc.mn]  1  . [c,mn] 

»U,U,L  P[a,U]  2  IZ,U1. 


+  H 


-[c,mn] 

s[a,kl]Ptc>”nl  ’ 


(23) 


Secondary  equation, 


’  [c,mn] 


-[c,mn]  r  \  *  -lc,mn] 

U*S]  “  j  ] 3 1 J  ,n]s[c,mn]  *  F[j,n] 


+  l  l 

a  [a,kl] 


-  H' 


:  .  5ic»mni 

n[a,kl]  lc*mnJ  ta’kl] 


+  H 


pI-.klJ  Pta,kl] 


+  H 


[a,kl] 


,  ..  ,  -[c,mn]  [c,mn] 

|Jrc,mn]  *  gs[a>kl]  2  [a.klj 


(24) 
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Normal  equation, 


#  r[c,mn]  c  c  *  rfc.mn] 

[c,nm]  [1,5]  j  Ij*n]n[c,mn]  *  F[j,n] 


+  I  I  -  h;  ,  •  g!c’”"! 

a  U.U)  n[..kl]  I'-""1  [a>kl) 
+  H  n,  ,  .  i[c'mnl 

p[..U]  Spia.U] 


-[c,mn] 

Ic,,ml  '  “V.klJ 


In  Eqs.  23,  24,  and  25,  the  following  definitions  have  been  used: 


[c,mn]  Used  to  indicate  that  the  test  point  Is  at  the  center  of 
patch  mn  on  region  c. 

[a,kl]  Evaluation  region,  indicates  that  the  integration  vari¬ 


[c,mn] 

[a,kl] 


able  x  runs  over  patch  kl  on  region  a. 

Kronecker  delta,  has  value  1  when  a  ■  c  and  m  -  k,  n  - 
1;  0  otherwise. 

Indicates  that  the  surface  integration  runs  over  field 
boundary  j . 
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fa* 


Gfa  kll  *  ♦"  do 

la»telJ  ta.kl] 


[c,mn]  _  £  p(V$  .  n)  da 

p[a,kl]  [a,kl] 


-U,mn]  _  r  s(7$  •  n)  da 

8[a,kl]  la»klJ 


We  now  have  as  many  equations  as  there  are  surface  patch  unknowns,  and 


1 

M 


The  testing  points  [c,mn]  in  Eq.  26  should  be  fairly  uniformly  distrib¬ 
uted  on  each  field  boundary  to  optimize  the  condition  number  of  the 
resulting  system  of  equations.  If  you  have  K  modes  on  each  port,  the 
number  of  test  points  L  on  each  port  required  to  determine  a  solution  is 
given  by  L  >  K/3;  L  must  be  an  integer.  The  factor  of  3  results  from 
the  fact  that  Eq.  26  is  a  vector  equation,  so  three  scalar  equations  are 
generated  for  each  test  point. 

The  equations  that  result  from  testing  on  the  conducting  bound¬ 
aries,  Eqs.  23,  24,  and  25,  and  those  due  to  testing  on  the  field  bound¬ 
aries,  Eq.  26,  form  an  overdetermined  linear  system  of  equations  and 
will  later  be  referred  to,  respectively,  as  type  A  and  type  B  equa¬ 
tions.  If  no  approximations  had  been  made  obtaining  these  equations, 
they  would  be  perfectly  consistent  and  a  unique  solution  satisfying  all 
the  equations  would  be  possible.  But  the  pulse  basis  expansions  are  not 
smooth  and,  hence,  cannot  perfectly  mimic  the  physically  exact  solu¬ 
tion.  As  a  result,  our  system  of  equations  is  inconsistent.  In  a 
strict  mathematical  sense,  the  system  of  equations  developed  here  has 
no  solution;  that  is,  there  exists  no  solution  vector  that  can  perfectly 
satisfy  every  equation  in  the  system.  However,  from  a  practical  point 


of  view,  the  system  has  minimum  residual  error  solutions  where  the 
residual  error  is  measured  by  some  suitable  vector  norm. 

Classically,  a  least-squares  residual  error  solution  can  be 
obtained  by  premultiplying  the  overdetermined  system,  Ax  *  b,  by  the 
conjugate  transpose  of  A,  A+,  and  solving  the  resulting  square  system  of 
equations  using  a  standard  method.  This  results  however  in  squaring  the 
condition  number  of  the  original  system  of  equations  and  should  be 
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avoided,  since  moment -net hod  matrices  are  notoriously  ill-conditioned  in 
the  first  place*  Methods  for  finding  the  pseudo-inverse  of  A  that  are 
based  on  finding  its  singular  value  decomposition,  which  is  nondependent 
on  premultipllcatlon  by  A+,  and  are  available  as  library  subroutines 


which  should  work.  But  they  run  at  substantial  overhead  cost  and 
require  a  significant  amount  of  additional  storage. 

For  the  results  in  Section  3,  a  variant  on  the  iterative  algorithm 
called  ART  was  employed.  The  synonym  ART  stands  for  ’’algebraic  recon¬ 
struction  technique’’  and  was  first  used  extensively  on  early  X-ray  CT 
scanners;  it  was  the  algebraic  tool  used  to  reconstruct  cross-section 
profiles  from  X-ray  attenuation  data,  and  hence  its  name.  This  algo¬ 
rithm  and  some  variants  thereon  are  explained  in  detail  in  Appendix  B. 
Briefly,  it  is  a  row  by  row  steepest  descent  algorithm  that  can  be 
relaxed  to  give  minimum  norm  error  solutions.  Being  a  "row  action” 
method,  ART  can  be  implemented  with  just  one  row  of  matrix  elements  in 
the  computer  core  at  a  time,  keeping  all  other  rows  on  disc  space. 
Because  disc  access  time  is  slow,  it  is  best  to  read  in  as  many  coeffi¬ 
cient  rows  as  possible,  so  that  the  number  of  disc  access  requests  is 
minimized.  Although  ART  exhibits  only  a  linear  rate  of  convergence, 
guarantee  of  convergence  does  not  depend  on  degree  of  diagonal  dominance 
as  with,  for  example,  the  SOR  algorithm;11  indeed,  ART  converges  almost 
unconditionally. 

If  the  system  of  equations  developed  in  this  section  were  per¬ 
fectly  consistent  and  very  well  conditioned,  there  would  be  no  need  to 
look  for  any  more  help  in  finding  the  solution  because  all  the  physics 
would  be  represented  and  mimicked  by  the  equation  system.  But,  as  was 
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noted  earlier,  com  physics  has  been  lost.  To  some  extent,  we  can 
reintroduce  the  physics  by  asking  our  approxiaate  solution  have  at  least 
sosm  property  In  common  with  the  exact  physical  solution.  This  takes 


the  fora  of  enforcing  a  conservation  of  power  constraint  on  the  solu¬ 
tion.  Mathematically  we  require  that  all  the  propagating  power,  real 
power  flow,  in  the  reflected  and  transaltted  aodes  oust  be  the  saae  as 
the  power  carried  by  the  incident  aode.  Since  the  normal  mode  functions 
in  Eq.  18  have  all  been  normalised  to  carry  the  saae  aaount  of  power,  we 
can  write. 


“I 


where  0  Is  the  suaaing  Index  over  only  the  propagating  aodes  on  port 

j.  This  constraint  was  enforced  at  the  end  of  each  ART  update  of  the 

entire  systsca  by  a  constant  phase  adjustaent.  The  update  takes  the 

(k) 

following  fora:  Let  ar.  ,  be  the  predicted  node  amplitudes  at  the  end 

11, h  J 

(k+) 

of  the  kth  ART  pass,  then  adjust  the  aode  aaplltudes  to  a..  ,  before 

1 » n  j 

starting  the  next  ART  pass  by  using 


(k+) 

•li.n] 


ly  y  a<k>  a<k>‘ 
|j  . . . 


V  7 

Doing  this,  the  values  of  a.  .  always  satisfy  Eq.  27  without  altering 

l  i  »h  J 

the  phase  of  the  initial  values.  Enforcing  this  constraint  also  tends 
to  accelerate  the  overall  convergence  rate  for  the  algoritha  by  reducing 
the  solution  space  to  only  that  class  of  coefficients  that  satisfy  Eq. 


-  24  - 


W 


*  >  V 
vv 
v> 
V«> 


.W 


■'.y.v, 

•\V 


A  A  «.•  .  A  ■  A  J*.  A  *  .  -  -  •  .  •  .  .  A  ■  o  .  A  A  -  ,  *  A"  .  A  A  AAA 

A  AAA  J  ^  **  V  1  *  '‘v'*  V  A  A*  A/'  -  *VaV,V.  •*.**‘,  ^A/V'v^  s 


777? 


27.  Other  constraints*'  uy  also  ba  used;  if,  for  example,  you  have 
symmetric  ports,  than  constraints  on  the  solution  based  on  that  symmetry 


should  be  used. 

As  an  aside,  Splelman13  used  a  foraulatlon  stellar  to  the  one 
described  here  (but  In  two  dimensions)  to  find  the  cutoff  resonances  of 
arbitarary  cross-section  waveguides.  When  there  are  no  ports,  Eqs.  23, 
24,  and  25  reduce  to  a  square  matrix  system  of  the  form  Ax  •  0;  this  is 
an  eigenvalue  problem  for  resonances  In  a  three-dimensional  cavity. 
Splelman  found  the  eigen  or  resonance  values  by  searching  for  frequen¬ 
cies  where  the  magnitude  of  the  determlnent  of  A  haa  mini  mums.  Comput¬ 
ing  a  matrix  determinant  Is  an  N  cubed  operation,  where  N  is  the  edge 
dimension  of  the  matrix.  Embedding  an  N  cubed  operation  into  a  search 
algorithm  means  slow  business.  This  is  also  a  nonlinear  eigenvalue 
problem.  Since  all  the  matrix  elements  vary  in  a  complex  way  with 
frequency,  it  cost  N  squared  numerical  integrations  over  surface  patches 
to  update  A  to  a  new  frequency.  Splelman  gets  away  with  all  this 
because  his  matrix  size  is  small,  only  21  by  21  elements.  It  seems 
unlikely  that  a  practical  algorithm  for  finding  eigen  frequencies  in 
three-dimensional  cavities,  based  on  this  approach,  would  be  possible. 
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IV.  RESULTS 


This  section  gives  nuaerical  results  for  s  few  staple  test  esses 


based  on  Eqs.  23.  24,  23,  and  26  In  Section  III.  The  job  was  broken 


down  into  three  sections.  A  Fortran  code  was  written  to  iapleaent 


each.  The  first  code  allows  one  to  build  s  data  file  specifying  arbi¬ 


trary  rectangles  in  three  diaenslons.  Discontinuities  that  could  be 


built  using  rectangles  were  the  only  type  considered.  Then,  for  each 


wavelength,  a  second  code  computed  on  the  order  of  N  squared  coupling 


coefficients,  where  N  is  the  total  nuaber  of  unknowns.  The  last  code 


'solved'  the  overdete rained  systea  of  equations  using  a  variant  on  the 


ART  algor it  ha.  These  Fortran  prograas  are  documented  in  Appendix  C. 


They  were  all  developed  on  the  Could  9080  coaputer  systea  at  the  Univer¬ 


sity  of  Utah  College  of  Engineering. 


Writing  a  code  capable  of  specifying  more  arbitrary  shapes  is 


certainly  possible  but,  froa  an  laplcsientation  point  of  view,  this  job 


of  specifying  and  organizing  the  surface  geoswtry  is  the  aost  diffi¬ 


cult.  Once  the  surface  shapes  are  specified  and  ordered  by  some 


accounting  systea,  the  calculus  follows  aechanically  froa  the  equations 


of  Section  III.  The  purpose  here  was  only  to  find  out  if,  and  then  how 


well,  a  formulation  of  this  type  could  be  expected  to  work.  So  while  a 


code  Halted  to  only  discontinuities  that  could  be  built  from  rectangles 


is  of  no  practical  value,  it  was  entirely  adequate  for  answering  many 


basic  questions  concerning  the  performance  of  the  algorithm. 


Data  are  presented  for  three  test  cases.  The  first  two  are  very 


simple,  a  straight  section  of  uniform  waveguide  as  a  two-port  junction 
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ud  a  straight  section  of  uni f on  waveguide  with  a  perfect  short  on  one 
side.  These  test  cases,  although  trivial,  are  iaportant  because  exact 
analytic  results  are  available  for  comparison.  They  were  used  Initially 
to  debug  the  Fortran  codes,  and  later  to  help  answer  such  questions  as 
how  many  surface  patches  are  needed  along  the  length  and  across  the 
width  of  the  waveguide  at  a  given  operating  wavelength  in  order  to 
retain  an  accurate  solution  without  using  an  excessive  nuaber  of 
patches.  For  the  «tralght  section  of  unifor*  waveguide,  results  are 
given  with  and  without  the  conservation  of  power  constraint  enforced. 

An  asymmetric  Iris  In  a  straight  section  of  uniform  waveguide  was 
considered  as  the  third  test  case.  These  results  were  compared  with 
analytic  results  found  using  approxiswte  susceptance  data  available  in 
the  Microwave  Engineers'  Handbook,14  Vol.  1,  p.  81.  For  the  reflected 
phaae  and  amplitude,  measured  data  were  also  obtained  using  a  slotted 
line  for  further  comparison.  Data  are  also  given  in  this  case  showing 
the  equation  consistency  (average  residual  error  per  equation)  as  a 
function  of  wavelength  for  the  solutions  presented. 


Case  1.  Straight  Shot  Data 


In  this  case,  the  TEjq  mode  In  a  rectangular  waveguide  is  numeri¬ 
cally  "propagated”  for  a  distance  of  one  tenth  of  the  free  space  cutoff 

wavelength,  X  ,  along  the  guide  axis.  In  Fig.  3,  the  cross-section 
c 

dimensions  were  a  -  0.5  X^  and  b  ■  0.2  *  ,  or  equivalently  by  0.5  by  0.2 
cutoff  units. 
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Fig.  3.  Waveguide  cross  section  for  straight  shot  data 


The  surface  patches  were  such  that  we  had  12  subdivisions  along  a, 
2  along  b,  and  2  sore  along  the  axis.  Each  field  boundary  was  tested  at 
24  uniformly  distributed  points.  This  resulted  in  168  surface  patch 
unknowns  and  2  field  boundary  unknowns.  Since  the  field  boundaries  were 
tested  so  many  times,  we  had  312  complex  equations. 

Two  sets  of  results  are  given.  Figs.  4  and  5  represent  the  solu¬ 
tion  obtained  without  using  the  conservation  of  power  constraint,  and 
Figs.  6  and  7  show  the  solution  with  the  constraint  enforced,  as  out¬ 
lined  in  Section  Ill .  On  all  the  plots,  the  abscissa  gives  the  operat¬ 
ing  wavelength  in  cutoff  wavelength  units.  The  upper  curve  in  Figs.  4 
and  6  represents  the  computed  transmission  amplitudes,  while  the  lower 
curves  on  these  plots  are  the  computed  reflection  amplitudes.  Clearly, 
the  exact  results  should  be  unity  transmission  and  zero  reflection 
amplitudes.  Figures  5  and  7  are  the  transmission  phase  with  respect  to 
the  input  plane  for  the  unconstrained  and  constrained  cases,  respect¬ 
ively.  On  these  plots,  the  marked  curves  are  from  the  computed  data, 
while  the  solid  unmarked  curves  represent  the  exact  phase  given  by 


where  the  operating  wavelength  is  giver,  by  X  ■  R  X  . 

c 
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Computed  transmission  and  reflection  amplitudes 
for  the  straight  shot  case  without  the  conserva 
tion  of  power  constraint  enforced. 
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Fig.  5.  Computed  versus  theoretically  exact  phase 
for  the  transmission  data  in  Fig.  4. 
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Fig.  6.  Computed  transmission  and  reflection  amplitudes 


for  the  straight  shot  case  with  the  conservation 
of  power  constraint  enforced. 


Fig.  7.  Computed  versus  theoretically  exact  phase 
for  the  transmission  data  in  Fig.  6. 
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From  these  data.  It  can  be  seen  that  the  conservation  constraint 
made  very  little  difference  in  the  solution.  The  transmission  phase  and 
amplitude  are  slightly  more  accurate  with  the  constraint  enforced,  but 
the  reflection  amplitude  was  a  bit  more  accurate  without  the  con¬ 
straint.  Note,  the  overall  accuracy  in  the  unconstrained  solution 
indicates  that  the  system  of  equations  used  has  retained  most  of  the 
physics  despite  the  expansion  function  approximations.  One  of  the  most 
remarkable  aspects  of  this  solution  was  the  ability  of  the  numerical 
algorithm  to  accurately  track  the  phase  all  the  way  from  the  low  fre¬ 
quency  end,  where  the  guide  wavelength  is  much  longer  than  the  free 


space  wavelength,  up  to  the  high  frequency  end,  where  the  guide  wave¬ 
length  and  the  free  space  wavelength  are  not  much  different. 


Case  2.  Straight  Shot  with  Perfect  Short 

For  this  case,  the  TE^q  mode  in  a  rectangular  waveguide  is  numeri¬ 
cally  reflected  off  a  perfect  shoring  plane  0.1  away  from  the  refer¬ 
ence  field  plane.  The  guide  cross  section  was  the  same  as  for  Case  1. 
Again,  12  subdivisions  were  used  along  a,  2  along  b,  and  2  along  the 
axis.  The  shorting  rectangle  was  subdivided  into  3  by  12  patches  along 
its  short  and  long  axes,  respectively.  The  field  boundary  was  tested  at 
36  uniformly  distributed  points.  This  results  in  276  surface  patch 
unknowns  and  1  field  boundary  unknown,  for  a  total  of  277  complex 
unknowns.  The  number  of  complex  equations  as  a  result  of  testing  is 


3 
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Again,  two  sets  of  results  are  given.  Figures  8  and  9  give  the 
amplitude  and  phase  for  the  unconstrained  solution,  and  Fig.  10  gives 
the  phase  for  the  constrained  solution.  In  the  case  of  the  constrained 
solution,  the  amplitude,  as  a  result  of  the  constraint,  is  exactly 
unity,  and  so  was  not  plotted.  Both  phase  plots  indicate  the  numeri¬ 
cally  predicted  solutions  with  marks;  the  unmarked  curves  are  the  exact 
phase  results  which  are  gxven  by 


♦  -  -  180°  -  72° 


1  -  R 


where,  as  before,  R  is  the  abscissa  value  on  these  plots.  The  abscissa 
value  is  related  to  the  physical  wavelength  by  X  »  R  X^. 

While  the  phase  tracking  is  not  bad  for  the  unconstrained  solu¬ 
tion,  the  reflection  amplitude  (Fig.  8)  deviates  away  from  unity  by  as 
much  as  10  percent  at  the  short  wavelength  end,  R  ■  0.4.  In  the  region 
where  the  data  overlap  the  phase  accuracy  of  the  constrained  solution  is 
a  bit  better  than  for  the  unconstrained  solution,  and  is  not  bad  all  the 
way  down  to  R  ■  0.2.  Of  course,  as  was  previously  mentioned,  the 
reflected  amplitude  is  perfectly  unity  for  the  constrained  solution. 
Although  the  unconstrained  solution  mimics  the  exact  solution  rather 
well,  the  value  of  the  constraint  condition  as  a  means  of  improving  the 
solution  is  much  clearer  here  than  in  Case  1. 


-  32  - 


g 


wmmmmmwmmmmmmzmMi 


1 

A 


■Transmission 


Wave  length 

Reflection  amplitude  for  the  straight  shot  case  with 
a  perfect  shorting  plane  on  one  side  and  also  without 
enforcing  the  conservation  of  power  constraint. 
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Fig.  9.  Computed  versus  theoretically  exact  phase 
for  the  reflection  data  in  Fig.  8. 
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Fig.  10.  Computed  versus  theoretically  exact  phase  for 

the  straight  shot  with  shorting  plane  case  with 
the  conservation  of  power  constraint  enforced. 


Case  3.  Straight  Shot  with  Asymmetric  Iris 


In  this  case,  an  asymmetric  iris  is  placed  0.1  away  from  the 
input  port,  and  the  output  port  is  0.1  on  the  other  side  of  the 
iris.  The  iris  is  assumed  to  be  perfectly  conducting  and  to  have  zero 
thickness  (see  Fig.  11).  The  waveguide  cross  section  is  the  same  as  for 
Case  1.  This  geometry  was  subdivided  using  16  divisions  along  a,  2 
along  b,  and  4  along  the  axis  from  port  to  port.  The  iris  was  subdi¬ 
vided  using  3  by  4  patches  for  6/a  *  0.75,  and  3  by  6  patches  for  6/a  * 
0.625.  Note  that  there  are  two  iris  rectangles,  one  for  each  side. 
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Fig.  11.  Geometry  for  asymmetric  iris. 


The  dominant  TEjq  mode  was  used  to  drive  the  system,  while  the 
first  8  TEjjq  modes  were  considered  in  the  field  expansions  at  each  field 
plane  or  port.  Each  field  boundary  was  tested  at  32  uniformly  distrib¬ 
uted  points.  Taken  all  together,  this  resulted  in  504  surface  patch 
unknowns  and  16  field  boundary  unknowns,  for  a  total  of  520  complex 
unknowns.  The  number  of  complex  equations  generated  by  testing  was 
696.  This  is  for  the  case  6/a  -  0.75.  In  the  case  6/a  -  0.625,  we  had 
556  complex  unknowns  and  732  complex  equations.  Computation  times  were 
on  the  order  of  25  minutes  per  wavelength  to  generate  the  coupling 
coefficients,  and  10  minutes  per  wavelength  to  solve  the  resulting 
system  of  equations.  These  times  are  for  the  Gould  9080  computer 
system  when  there  are  no  other  users  on  board. 

The  conservation  of  power  constraint  was  enforced  for  all  of  the 
computed  curves  for  Case  3.  The  numerically  computed  solutions  are 
marked  by  ”o"  for  TE^q  results,  and  by  "*”  for  results.  Reflection 
phase  measurements  were  made  using  a  slotted  line,  and  these  results  are 
denoted  using  an  "I"  symbol  on  the  reflection  plots.  The  solid  unmarked 
curves  are  results  predicted  using  published  susceptance  data,  they  are 
for  TE,n  comparison. 
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Maracuvitz15  indicates  that  susceptance  data  of  the  type  used 
here  have  been  derived  by  the  “equivalent  static  method,”  which  uses  the 
static  aperture  field  for  the  two  lowest  modes.  In  terms  of  these  data, 


I 


the  complex  reflection  amplitude  is  given  by 

a  ^i8L 

and  the  complex  transmission  amplitude  by  i 

+  e~2jBL 
a  (2L)  “  2(Yq/BJ  -  j 

where  the  ratio  Yq/b  is  found  on  page  81  of  the  Microwave  Engineers' 

Handbook, and  8  «  2*/A  .  For  all  the  data  in  this  section,  L  * 

g 

0.1  *  .In  the  limiting  case  of  a  perfect  short,  YQ/B  goes  to  0,  and  the 

-21BL 

equation  for  the  reflection  amplitude  gives  a  (o)  -  e  ,  which  is  the 

correct  amplitude,  but  the  phase  is  wrong  by  180°.  One  can  simply 
multiply  the  equation  by  -1  to  correct  the  deficiency,  but  rechecking 
the  derivation  did  not  warrant  this.  The  reflection  phase  curves  were 
all  plotted  with  this  180°  discrepancy  removed.  Shifting  these  curves 
180°  brings  them  into  closer  agreement  with  both  the  measured  and  numer¬ 
ically  predicted  phase  curves. 

Figures  12-15  are  the  computed  solutions  for  6/a  «  0.75.  In  this 
case,  agreement  between  the  numerical  solution  and  the  solution  found 
using  susceptance  data  is  really  quite  good,  especially  for  the  trans¬ 
mission  data.  Figures  16-19  give  solutions  for  6/a  -  0.625.  Here  the 
agreement  is  less  striking.  Keep  in  mind  that  the  susceptance  data  are 


a  '■o;  “  2T 
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not  exact.  Compare  the  reflection  amplitude  curves  in  Figs.  12  and 
16.  Note  that  the  susceptance  data  solutions  have  changed  quite 
markedly,  while  the  numerically  predicted  solutions  differ  more 
modestly.  The  numerically  predicted  solutions  also  appear  to  vary  in  a 
more  realistic  way  near  the  transition  region,  where  R  ■  1/2.  Although 
not  conclusive,  the  measured  reflection  phase  angle  is  closer  to  the 
numerically  predicted  curves  than  to  the  susceptance  data  curves,  even 
when  these  curves  are  corrected  by  180°,  as  previously  mentioned. 

Figures  20  and  21  give  a  measure  of  the  consistency  in  the  equa¬ 
tion  set  at  that  point  in  the  Iterative  process  where  the  solution  was 
taken.  The  ordinate  on  these  graphs  gives  the  average  error  per  equa¬ 
tion.  The  reason  for  this  inconsistency  was  detailed  in  Section  III;  it 
is  primarily  due  to  the  pulse  basis  approximation.  These  graphs  show 


that  the  equation  sets  are  more  consistent  at  the  high  frequency  end 
than  at  the  low  frequency  end.  Although  this  was  not  expected  a  priori, 
one  postulation  is  that  at  the  low  frequency  end,  the  guide  wavelength 
is  much  different  from  the  free  space  wavelength,  and  the  effect  of  the 
waveguide  walls  is  much  stronger  making  the  surface  patch  approximation 
more  important.  On  the  other  hand,  at  the  high  frequency  end,  the  guide 
wavelength  and  the  free  space  wavelength  are  nearly  the  same,  so  the 
effect  of  the  guiding  walls  is  less  critical. 
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Fig.  20.  Average  error  per  equation,  case  6/a  -  C 


Wavalangth;  Seala  “  i.E- 
Fig.  21.  Average  error  per  equation,  case  6/a  -  0 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  paper  has  presented  a  three-dimensional  multiport  formulation 
for  waveguide  scattering.  Simple  test  computations  were  made  to  verify 
the  validity  of  the  formulation.  Test  data  were  found  to  be  in  good 
agreement  with  results  obtained  using  other  methods  for  TEjq  compari¬ 
son.  In  the  case  of  the  asymmetric  iris,  results  were  obtained  for 
higher  order  TE  modes  as  well;  these  results  are.  of  dubious  accuracy,  as 
no  other  data  were  available  for  comparison.  It  is  suspected  that  the 
higher  mode  accuracy  will  diminish  rapidly  because  the  patch  size 
becomes  comparable  to  the  rate  of  variation  in  these  modes. 

Although  the  formulation  presented  here  is  quite  generally  valid, 
it  may  be  more  expensive,  both  in  terms  of  storage  cost  and  in  terms  of 
computational  intensity,  than  a  corresponding  three-dimensional  finite 
element  approach.  Computationally,  one  must  compute  on  the  order  of  N 
squared  numerical  integrations  at  each  new  wavelength.  Because  the 
solution  algorithm  is  iterative,  the  same  coefficients  are  needed 
several  times,  and  hence  all  N  squared  complex  numbers  must  be  saved. 
The  recomputation  cost  is  simply  too  high  to  do  otherwise.  This  repre¬ 
sents  a  storage  cost  proportional  to  the  surface  area  of  the  enclosed 
volume  squared,  or  equivalently  proportional  to  the  volume  to  the  four- 
thirds  power.  Finite  element  methods  would  require  storage  proportional 
to  the  volume.  So  a  formulation  of  this  type  may  not  be  competitive 
with  the  finite  element  method  for  internal  scattering  problems.  The 
formulation  of  this  paper  may,  however,  be  valuable  for  determining  the 
scattering  characteristics  of  slots  in  waveguides;  that  is,  junctions 


that  involve  radiation  into  free  space.  Here,  boundary  conditions 
external  to  the  waveguide  must  also  be  taken  into  account.  Unbounded 


problems  pose  more  difficulties  for  finite  element  formulations  than  for 
surface  integral  equation  approaches. 

The  formulation  in  this  paper  could  be  pushed  into  a  category  of 
useful  mathematical  software  if  some  of  the  following  problems  could  be 
solved: 

1.  A  generalized  surface  generation  routine  could  be  developed 
that  would  allow  the  engineer  to  rapidly  define  new  surface 
geometries.  This  might  be  done  using  a  triangular  net  bound¬ 
ary.  Each  net  point  could  then  be  moved  about  in  three  dimen¬ 
sions  to  form  a  wide  variety  of  surface  shapes. 

2.  Because  the  Integral  operator  is  a  “smoothing"  operator,  we 
got  away  with  using  a  pulse  basis  expansion  on  the  conducting 
boundaries  for  the  simple  test  cases  of  this  paper.  However, 
it  is  suggested  for  more  complex  geometries  and  increased 
solution  accuracy  that  a  smooth  C°  basis  set  be  used.  For  the 
triangular  net  above,  the  rooftop  functions  described  in 
Section  III  would  form  an  appropriate  basis  set. 

3.  For  the  asymmetric  iris  data,  the  bulk  of  the  computation  time 
was  tied  up  in  computing  the  coupling  coefficients  and  not  in 
solving  the  resulting  system  of  equations.  Algorithms  capable 
of  rapidly  and  accurately  approximating  these  coefficients 
would  be  valuable. 

4.  Many  important  problems  involve  at  least  one  plane  of  sym¬ 
metry;  if  this  symmetry  could  be  exploited,  the  number  of 
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IJ 


APPENDIX  A 


NORMALIZED  MODE  FUNCTIONS  AND  EXPANSIONS 


Part  1 


General  relations  for  EM  waves  on  a  uniform  cylindrical  waveguide- 
ing  system  with  care  taken  to  preserve  sign  changes  as  a  result  of 
assumed  propagation  direction. 

In  a  uniform  waveguiding  system,  that  is,  a  system  having  no 
physical  or  material  variation  along  the  z  axis,  it  turns  out  that  we 
have  simplified  relations  for  the  transverse  components  of  any  mode 

A 

supported  by  the  structure  in  terms  only  of  the  axial,  z  directed  field 
components. 

To  this  end,  expand  the  field  vectors  and  the  del  operator  into 
their  transverse  and  axial  components. 


A  5  A  +  A 
t  z 


7  i  7  +7 

t  z 


With  these  definitions,  note  that 


7  x  A  *  7  x  A.  +  7  x  A  +  7  x  A„ 
t  t  t  z  z  t 


Using  this  identity  in  Maxwell’s  equations,  Eqs.  10  and  11,  and  iden¬ 
tifying  component  directions,  we  obtain  two  relations  for  the  transverse 
components  of  the  equations. 
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^  v  V  V  '>  */  VV  V  V  *  a  V  •  •  *•*■'•*  *  •  *  i 


■..VVV’V 
■\W  - 


v  x  h  +7  XH-  lueE 
t  z  z  t  t 


(A. 2) 


Since  the  structure  is  uniform  along  z,  we  can  assume  that  mode  m  will 
propagate  with  spatial  variation  that  looks  like 


z 

A(x,y,z)  =  A(x,y)  e 


where  the  -  sign  is  for  propagation  in  the  +  z  direction  and  the  +  sign 

A 

for  propagation  along  the  -  z  direction.  With  this  understanding,  note 


V  x  A  -  ?  jfi  (z  x  a  ) 
z  mt  nr 


Using  this  result  in  Eqs.  A.l  and  A. 2, 


-  V  x  e  +  18  (z  x  I  ) 
in  t  m  J  nr  m  J 


(A. 3) 


IojeE  «  V_  x  H  +  IS  fz  x  H  1 
J  m  t  m  J  nr  m  J 
t  z  t 


(A. 4) 


Substituting  Eq.  A. 4  into  Eq.  A. 3  to  eliminate  the  transverse  component 
of  E  and  using  the  identities 


r,  -r.  *.r.  •>.  -v*.  ■*. 

/,v,  /.Vv.v  .'.v, 


/  /  /  .  t*  .* 


LWWMWf  HMIUMMA  U 


A 

z  X  (v  X  A  )  ■  7  A 
v  t  ZJ  t  z 

(z  X  At)  -  -  At 


Z  X 


we  obtain 


2  2  —  — 

(k  -  B  )H  =  joieV  x  e  +  jB  7  H 
v  b,  t  in  mtm 

t  z  z 


(A. 5) 


Similarly,  we  could  eliminate  the  transverse  components  of  H  by  substi¬ 
tuting  Eq.  A. 3  into  Eq.  A. 4  to  obtain 


(k2  -  02)I 

'  ra'  in 


-  jup7^  x  H  +  jB  7  E 
J  t  m  J  m  t  m 


(A. 6) 


2  2 

Here,  k  2  u>  pe.  Equations  A. 5  and  A. 6  are  standard  results  except  that 


the  +  signs  are  not  often  Included;  usually,  just  the  -  sign  is  used  for 

A 

propagation  in  the  +  z  direction.  This  distinction  is  important  for 
mode  expansions  involving  propagation  in  both  directions. 


In  the  case  of  TE  modes,  Em  =  0  and  Eqs.  A. 5  and  A. 6  reduce  to 

z 


H 


m 


7^  H 


„2  t  m 
t  k  -  B  z 

m 


-jn»P 


mt  k2  -  B2 


ra 


[7<  v  ;J 


(A. 7) 


For  TM  modes,  Hm  “  0  and  Eqs.  A. 5  and  A. 6  reduce  to 
z 


-  48  - 


& 
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-  M  " 

*  j»  v  V, 


X^.\VN\VV/ V  \  ,  •  ►  W/v  *// 

. , 


•v*.; 

mm 


These  equations  show  that  the  axial  field  component  in  a  cylindri¬ 
cal  waveguide  plays  a  role  similar  to  that  of  a  scalar  potential  in 
general  EM  field  theory.  The  fields  in  a  cylindrical  system  are  com¬ 
pletely  determined  once  the  axial  fields  are  known. 

Since  the  electric  and  magnetic  fields  are  solutions  to  vector 
wave  equations,  Eqs.  14  and  15,  the  axial  fields  must  he  solutions  to 


( V2  +  y2) 
t  mJ 


(A. 9) 


2  2  2 

where  y  =  k  -  S  •  Solutions  to  Eq.  A. 9  are  subject  to  the  following 
m  m 

boundary  conditions: 


n  •  VH  -  0  for  TE 
m 

z 


E  *  0  for  TM 
m 

z 


where  n  is  a  unit  normal  to  the  conducting  surface  of  the  cylindrical 
waveguide. 
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Part  2 

Normalized  mode  functions.  In  Part  1,  the  basic  relatio)  ships 
between  field  components  in  a  cylindrical  system  were  outlined.  The 
electric  and  magnetic  mode  vectors  were  denoted  by  capital  letters  E  and 
H,  respectively.  In  this  section,  a  normalization  condition  is  given. 
The  normalized  mode  -vectors  will  be  denoted  using  lower  case  e  and  h. 
The  normalization  condition  used  here  is  that  each  mode  carry  1  watt 
real  power  above  cutoff,  and  1  watt  reactive  (imaginary)  power  below 
cutoff.  Explicitly,  integrate  the  Poynting  vector  over  the  guide  cross 
section  S  and  enforce  the  following  condition: 


f 


1 

[W] 

for 

both  TE  and  TM 

modes  above 

cutoff 

1  r  * 

2  J  em 

•  ds  *■  J 

m 

j 

tw] 

for 

TE  modes  below 

cutoff 

(A. 10) 

S  t 

t 

-j 

[WJ 

for 

TM  modes  below 

cutoff 

The  change  of  sign  on  the  reactive  power  below  cutoff  reflects  the  fact 


that  nonpropagating  TE  modes  appear 
modes  appear  capacitive. 

Using  the  relations  from  Part 
A.  10  can  be  enforced  in  terms  of  the 
use  Eqs.  A. 7  to  obtain 


* 


m 


For  TM  modes,  use  Eqs.  A. 8  yielding 


inductive,  while  nonpropagating  TM 

1,  the  condition  indicated  by  Eq. 
axial  fields  only.  For  TE  modes, 

1  [W]  above  cutoff 

(A. 11) 

1  [W]  below  cutoff 


m 
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1  [W]  above  cutoff 
-j  [W]  below  cutoff 


(A. 12) 


Once  the  axial  modes  have  been  normalized,  the  corresponding  transverse 
field  components  may  be  found  using  Eqs.  A. 7  and  A. 8. 

Arbitrary  fields  in  a  cylindrical  waveguide  can  now  be  written  in 
terms  of  the  normalized  mode  functions, 


H  ■  £  (ah  +  b  h  ) 

**  v  m  m  m  nr 
m 

where  m  ranges  over  all  modes,  both  TE  and  TM,  above  and  below  cutoff. 

A  A 

The  forward-backward  arrows  indicate  propagation  in  the  +  z  and  -  z 
directions,  respectively. 

Part  3 

Explicit  mode  functions  for  rectangular  waveguides.  In  a  rectan¬ 
gular  waveguide  with  cross-section  dimensions  a  and  b,  the  mode  func¬ 
tions  that  satisfy  Eq.  A. 9  with  the  appropriate  boundary  condition  for 
TE  modes  is 


with 


h 

z 


mn 


A 

mn 


xmx 

cos  -  cos 

a 


....;  but  not  both  zero 


m,n  *  0,  1,  2, 


Using  Eq.  A. 12  for  normalization. 


A  - 

mn 


2 

2vy  n 
mn  mn 


1 1/2 


kg*  ab 
mn 


If 


where 


1  [W]  above  cutoff 


j  [W]  below  cutoff 


n  * 

mn 


2  if  m  or  n  is  0 


4  otherwise 


For  TM  modes,  the  mode  function  that  satisfies  Eq.  A. 9  with  >  0  on 
the  boundary  is 


i  -  A  .In  IS  sin  1ST 
z  mn  a  b 

mn 


with 


2  /  mn  y  /  nn  \2 

^mn  "  \  a"~ )  +(-j 


m,n  ■  1,  2,  .... 


Using  Eq.  A. 13  for  normalization. 


A 

mn 


8vy2  j — 

mn  I  w 

kB  ab  |  r 
mn 


1/2 
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1  [ W]  above  cutoff 


v 


-j  fW]  below  cutoff 


APPENDIX  B 


COMPLEX  ART 

Consider  a  linear  system  of  equations. 


Ax  ■  b 


00 

Let  %  be  the  kth  approximation  to  the  solution  vector  x  and  note  that 
the  error  or  residual  in  the  ith  equation  in  the  system  can  be  written 
as 


where  a^  is  the  ith  row  vector  from  the  matrix  A,  and  b^^  the  ith  element 
in  the  column  vector  b.  The  gradient  of  the  residual,  steepest  descent 

direction,  can  be  computed  with  respect  to  a  change  in  the  approximation 

„  00 

vector  y  , 


This  says  that  the  gradient  of  the  residual  for  the  ith  equation  in  the 
system  is  the  transpose  of  the  ith  row  vector  in  the  matrix  A.  Now, 
adjust  the  approximate  solution  along  the  steepest  descent  direction 
using  the  update  form 


^(k+l) 


00  T 

z  ak  Si 
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The  scalar  is  chosen  to  make  r^  zero.  This  gives  the  condition 


0  ‘  **  [*(k>  ‘  5‘]  ‘  b‘ 


Solving  for  a^,  we  have 


T 

-i  -i 


The  basic  ART  update  for  real  systems  of  equations  is  then 


(k+1)  (k) 

y  -  y 


r(k)  aT 
ri _ -i 


-i  Si 


(B.l) 


Note  that  ART  adjusts  every  element  in  the  solution  vector  for  each 
equation.  When  Eq.  B.l  has  been  enforced  for  every  equation  in  the 
system,  a  complete  update  has  been  accomplished.  Usually,  several 
updates  are  required,  but  if  the  system  has  a  unique  solution,  then 


lim  -  x)  -  0 


Any  n  by  m  complex  system  of  equations  can  be  written  as  a  2n  by  2m  real 
system  to  which  Eq.  B.l  can  be  applied  directly.  Alternatively,  update 
formulae  equivalent  to  Eq .  B.l  can  be  derived  which  apply  directly  to 
complex  systems.  For  A  x  *  b  a  complex  system,  decompose  the  complex 
residual  into  its  real  and  imaginary  parts, 


ri  ’  °i  +  J6i 


a, B  real  valued 


SKy 

$ 

;.v;v 


if 

1W.& 


d 


i 

m 


i 
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Then  the  update  that  makes  zero  is 


(k)  + 

fk+1)  fk)  °i  -i 
Z  -  Z - “ 


(B.2) 


The  update  making  8^  zero  is 


(k+1) 


-i  -i 


.00  T 

.(k)  8i  -i 


-i  Si 


(B.3) 


where  the  following  notation  applies: 


a^  :  the  transpose  of  a^ 


a^  :  the  conjugate  transpose  of  a^ 


Equations  B.2  and  B.3  are  applied  by  computing  the  real  part  of  the 
residual  for  equation  i,  a^,  then  enforcing  Eq.  B.2.  Next  compute  the 
imaginary  part  of  the  residual  for  equation  i,  6^,  and  enforce  Eq .  B.3. 

A  variant  on  this  algorithm  that  I  call  new  ART  or  NART  applies  to 
systems  where  the  diagonal  element  is  larger  in  magnitude  than  other 
elements  in  the  same  row.  This  update  changes  only  the  ith  component 
of  jr  when  applied  to  the  ith  equation,  rather  than  all  the  elements 
in  ][*  The  NART  update  formulas  corresponding  to  Eqs.  B.2  and  B.3  are 
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where  the  elements  of  A  are  [a^J.  Equations  B.4  no  longer  reduce  the 
ith  residual  to  zero  at  each  update,  but  in  many  cases  this  update 
scheme  converges  in  fewer  iterations  and  it  costs  less  per  iteration. 
MART  is  the  update  scheme  used  by  subroutines  ARTA  and  ARTB  in  this 
report . 
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APPENDIX  C 


FORTRAN  FOREVER 

This  appendix  gives  documentation  for  the  three  Fortran  codes  used 
to  test  the  formulation  presented  in  Section  III.  These  codes  are  of  no 
Immediate  value,  but  they  are  Included  here  to  demonstrate  work,  success¬ 
fully  completed  and  also  because  they  may  be  useful  for  showing  how  such 
computations  may  be  organized. 

Program  rspec.f 

This  program  allows  for  the  specification  of  arbitarary  rectangles 
in  three  dimensions.  First,  the  plane  of  the  rectangle  is  specified  by 
a  point  in  the  plane,  Xj  ,  and  by  the  direction  of  the  "outward"  normal 

A 

to  the  plane,  n.  Here,  the  term  outward  means  outward  with  respect  to 
the  volume  that  is  enclosed  by  the  rectangles.  The  equation  of  the 

^  A  A 

plane  is  given  by  all  points  x  such  that  x  •  n  **  x^  •  n,  or  more  explic¬ 
itly  by 


n^x  +  n2y  +  n^z  »  n  •  Xj  (C.l) 

A  A  A  A 

where  n  ■  n^x  +  n2y  +  n^z«  The  second  vertex  of  the  rectangle  X2  can 
now  be  specified  using  only  two  of  its  three  components.  Since  specify¬ 
ing  two  of  the  three  variables  in  Eq.  C.l  completely  determines  the 
third  variable,  we  can  write  the  triplet  of  points  as: 
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'WVv 'vS. Vi .*.V*»*. ‘•Y -.-.y -V- , 


[(n  •  Xj  -  n2^2  ”  n3*2^nl*  ^2*  Z2^  ^or  ni  ^  n?>  n 


2  '  «2 9  ^ 


* 

[x^  (n  •  Xj  -  n^x2  ”  n3z2^n2'  z2^  for  n2  ^  ni  ’  n,a 


[x2,  y2,  (n  •  Xj  -11^2  -  n2y2)/n3]  for  n3  >  iij ,  n. 


(C.2) 


The  three  forms  above  are  entirely  equivalent  except  when  n^ *  n2»  or  n^ 
are  zero.  To  avoid  division  by  zero,  use  the  appropriate  version  of  Eq. 

A  A 

C.2;  note  that  at  least  one  component  of  n  is  not  0,  since  |nj  *  1.  The 
surface  primary  direction  can  now  be  defined  using 


p  ■  (*i  -  xx)x  +  (y2  “  yA)y  +  <z2  ~  zi): 


(C.3) 


The  third  vertex  x3  can  now  be  specified  by  entering  just  one  of  its 
three  components,  since  it  must  lie  in  the  plane  and  on  a  line  perpen¬ 
dicular  to  p,  i.e.,  s  •  p  *  0,  where  s  is  the  secondary  direction  given 


—  ^ 
8  "  (x3  -  Xj)x  +  (y3  -  yx)y  +  (z3  -  Zj)z 


(C.4) 


Enforcing  this  condition,  we  have  for  each  of  the  three  cases  in  Eq . 
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If  |a|  >  jbj,  use  X3  ■  (c/a)  -  (b/a^^.  enter  z^,  evaluate  x^,  then  use 
Eq.  C.l  to  find  y^,  otherwise  use  Z3  *  (c/b)  -  (a/b)x3,  enter  X3,  evalu¬ 
ate  Z3,  then  use  Eq.  C.l  to  rind  73. 

Case  3,  n^  >  n^,  n2 

ax^  +  by^  ”  c 

where 

3  "  x2  “  X1  ~  nl^Z2  "  Zl^/n3 

b  "  y2  "  yl  “  n2^Z2  "  Zl^/n3 

C  -  x,Cx2  -  X,)  +  y;(y2  -  y,}  ♦[*,-"•  -  x,) 

If  |  a  |  >  j  b |  ,  use  X3  -  (c/a)  -  (b/a)y3,  enter  y3,  evaluate  x3,  then  use 

Eq.  C.l  to  find  Z3 ,  otherwise  use  y^  *  (c/b)  -  (a/b)x^,  enter  x^,  evalu¬ 
ate  y^,  then  use  Eq.  C.l  to  find  z^. 

This  method  of  rectangle  specification  requires  a  minimum  amount 
of  input  from  the  user,  thereby  avoiding  specification  errors. 

Note  that  the  rectangle  vertex  information  should  be  ordered  so 

AAA  A  A 

that  p  x  s  -  n  ,  where  p  and  s  are  unit  vectors  in  the  directions  p  and 
s  from  Eqs.  C.3  and  C.4.  Any  point  on  the  rectangular  surface  can  now 
be  described  by  the  local  coordinates  (p,s).  The  transformation  from 
the  local  system  to  the  absolute  point  x  is  accomplished  by 


for 


(C.5) 


x  ■  Xj  +  pp  +  ss 


where 


P€[0,  Dp] 
s^[0,  Dg] 


D 

P 


(C.6) 


For  the  final  output  record  rspec.f  stores  the  following  informa¬ 
tion  about  each  rectangle: 

AAA 

Xj,  n,  p,  s,  Dp,  Np,  Ng 


where  D  and  D  are  the  primary  and  secondary  edge  lengths  from  Eqs. 
p  s 

C.6,  and  N  and  N  are  the  number  of  primary  and  secondary  subdivisions 
P  s 


along  each  edge 


PROGRAM  RSPEC 


This  program  allows  you  Co  specify  arbitrary  rectangles  in  three 
dimensions.  First  the  plane  of  the  rectangle  is  defined  by  the 
"outward"  normal  direction  and  the  coordinates  of  the  first  vertex, 
VI.  Next  the  program  will  ask  you  to  give  the  two  coordinates  need¬ 
ed  to  specify  the  second  vertex,  V2 ,  on  the  defined  plane.  Finally 
you  will  be  ask  to  enter  the  single  remaning  coordinate  required 
to  specify  the  third  vertex,  V3,  along  the  normal  to  V1-V2,  through 
VI,  and  on  the  defined  plane.  From  this  data  the  primary  and  second¬ 
ary  unit  directions  are  computed  for  later  use  by  surface  integra¬ 
tion  routines. 

REAL  VI (10, 3) ,UN( 10,3) ,UP( 10,3) ,US( 10,3) ,DP( 10) ,DS( 10) 

INTEGER  NP( 10 ) , NS ( 10 ) 

WRITEC6,*)  'Enter  the  number  of  rectangles  to  be  specified.' 

READ ( 5 , * )  NR 

WRITE(6,*)  ‘Enter  the  number  of  conducting  rectangles  out  of  total. 
READ ( 5 , * )  NRC 

WRITE ( 6 ,* )  ‘Specify  the  conducting  rectangles  first.' 

DO  100  I • 1 , NR 
WRITE ( 6 , * )  '  ' 

WRITEC 6 ,* )  'Data  entery  for  rectangle  number  *,I,'.' 

WRITEl 6 , * )  'Enter  normal  direction,  Nx,Ny,Nz.’ 

READC  5 , * )  X , Y , Z 
D-SQRT(X**2+Y**2+Z**2) 

UN(I ,1)-X/D 
UN (1,2) «Y/D 
UN (1,3) *Z/D 

WRITE( 6 , * )  'Enter  first  vertex,  X1,Y1,Z1.' 

READ ( 5 ,* )  XI ,Y1 ,Z1 

D-X1*UN(I,1)+Y1*UN(I,2)+Z1*UN(I,3) 

Test  for  free  variables. 

XN-ABS(UN(I,1) ) 

YN" ABS ( UN ( I , 2 ) ) 

ZN" ABS ( UN ( I , 3 ) ) 

IF( (XN.GE.YN) .AND. (XN.GE.ZN))  GOTO  60 
IF( (YN.GE.XN) .AND. (YN.GE.ZN) )  GOTO  30 

Case  1.  N3  >  N1  and  N2  . 

WRITE ( 6 , * )  'Enter  X2,Y2.' 

READC 5  ,  * )  X2,Y2 

Z2» ( D-UN (1,1) *X2— UN( 1,2) *Y2 )/UN(I,3) 

UP(I,1)-X2-X1 

UP(I,2)-Y2-Y1 

UP(I,3)-Z2-Z1 

A»UP( 1,1 )— UN ( I,1)*UP(I,3)/UN(I,3) 

B*UP( I ,2 )-UN( I , 2 )*UP( I ,3) /UN (I ,3) 

C"X1*UP( 1,1 )+Yl*UP( I , 2 )+( Zl-D/UN (1,3) )*UP(I,3) 

IF ( ABS ( A ) . GT . ABS ( B )  )  GOTO  10 

WRITEC6,*)  'Enter  X3. ' 

READC  5 , * )  X3 
Y3-C/B-(A/B)*X3 
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GOTO  20 


WRITE ( 6 ,  * )  'Enter  Y3.‘ 

READ(S,*)  Y3 
X3«C/A-(B/A)*Y3 

Z3-(D-UN(I,1)*X3-UN(I,2)*Y3)/UN(I,3) 

GOTO  90 

Case  2.  N2  >  N1  and  N3  . 

WRITE (  6  ,  *  )  'Enter  X2,Z2.' 

READ( 5 ,* )  X2,Z2 

Y2 - ( D-UN (1,1) *X2-UN ( I , 3 ) *Z2 ) /UN (1,2) 

UP (1,1) »X2-X1 
UP(I ,2)«Y2-Y1 
UP( I ,3  )»Z2-Z1 

A*UP(I,1) -UN ( I , 1 ) *UP( I , 2 ) /UN ( I ,2) 

B*UP( 1,3 )— UN (1,3) *UP (I,2)/UN(I,2) 

C-X1*UP( I , 1 )+Zl*UP( I , 3 )♦( Yi-D/UN( 1,2) )*UP(I,2) 
IF ( ABS ( A ) . GT . ABS ( B ) )  GOTO  40 

WRITE(6,*)  'Enter  X3.' 

READ( 5 , * )  X3 
Z3*C/B- ( A/B ) *X3 
GOTO  50 

WRITE ( 6 , * )  'Enter  Z3.' 

READ( 5 , * )  Z3 
X3**C/A-  (  B/A )  *Z3 

Y3«(D-UN(I,1)*X3-UN(I,3)*Z3)/UN(I,2) 

GOTO  90 

Case  3.  N1  >  N2  and  N3  . 

WRITE ( 6 , * )  'Enter  Y2,Z2. ' 

READ(5,*)  Y2,Z2 

X2« ( D-UN (1,2) *Y2-UN (I,3)*Z2)/UN(I,1) 

UP (1,1) "X2-X1 
UP( 1,2) "Y2-Y1 
UP(I ,3)-Z2-Zl 

A-UP( I , 2 )-UN( I ,2)*UP( I ,1)/UN( I ,1 ) 
B-UP(I,3)-UN(I,3)*UP(I,1)/UN(I,1) 

C»Y1*UP(I,2) +Z1*UP( I , 3 ) + ( Xl-D/UN (1,1) )*UP(I ,1 ) 
IF( ABS(A) .GT.ABS(B) )  GOTO  70 

WRITE( 6 , * )  'Enter  Y3 . ' 

READ( 5 , * )  Y3 
Z3-C/B-(A/B)*Y3 
GOTO  80 

WRITE( 6 , * )  'Enter  Z3. ' 

READ( 5 , * )  Z  3 
Y3-C/A-(B/A)*Z3 

X3« ( D-UN (I ,2)*Y3-UN(I ,3)*Z3)/UN( 1,1) 
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C  Compute  unit  directions  and  distances. 

90  DP( I ) “SQRTC  UP( I,1)**2+UP(I,2)**2+UP(I,3)**2) 
UP( I f 1 ) "UP( 1,11 /DP( I ) 

UP( I , 2 ) “UP( 1,2) /DP( Z ) 

UP( 1,3) »UP( 1,3 )/DP{ Z ) 

US(I,1)-X3-X1 
US(I ,2)-Y3-Yl 
US(I,3)-Z3-Z1 

DS( X ) »SQRT( US( 1,1 ) **2+US (1,2 ) **2+US( 1 , 3 ) **2 ) 
US ( I,1)»US(I,1) /DS( X ) 

US (1,2) »US (1,2) /DS ( I ) 

US(I,3)-US(I,3)/DS(I) 

VI ( I , 1 ) »X1 
Vl(I,2)-Yt 
V1(I,3)-Z1 


100 


110 


120 


WRITE(6,*)  'Enter  the  number  of  primary  subdivisions,  NP. ' 
READ! 5,*)  NP( I ) 

WRITE (  6  ,  *  )  'Enter  the  number  of  secondary  subdivisions,  NS. 
READ(5,*)  NS ( I ) 


Error  alowance,  in  case  of  incorrect  data  entry. 

WRITE ( 6 , * )  'Xf  data  entry  was  incorrect  type  *0“,  otherwise 

READ (5,*)  NANS 

IF ( NANS . EQ . 0 )  GOTO  5 

CONTINUE 


Store  point  data. 

OPEN ( 10 , FILE* ' rec . dat ' , STATUS- ' unknown ' ) 
REWIND ( 10 ) 

WRITEC 10 ,* )  NR,NRC 
DO  110  X-1,NR 

WRITE( 10 , * )  (V1(I,J) ,J-1,3) 

(UN( I , J) , J-l , 3 ) 

(UP(I,J) ,J-1,3) 
(US(I,J),J«1,3) 

DP( I ) ,DS( X ) ,NP( I ) , NSC  I ) 


WRITEdO,*) 
WRITEdO,*) 
WRITEdO,*) 
WRITE! 10,*) 
CONTINUE 
CLOSE (10) 


Test,  compute  vertices  from  data  as  a  check 
DO  120  I ■ 1 , NR 
WRITE!  6  ,  *  )  I 

WRITE( 6 , * )  (VI ( I , J) , J*1 , 3 ) 

X2-V1 ( I , 1 )+DP( I ) *UP( 1,1) 

Y2«V1 (1,2) +DP( X)*UP(I,2) 

Z2-V1 ( I , 3 )+DP( I ) *UP( 1,3) 

WRITE( 6 , * )  X2,Y2,Z2 
X3*V1 ( I , 1 ) +DS ( I ) *US ( I , 1 ) 
Y3-V1(I,2)+DS(I)*US(I ,2) 

Z3*V1 (1,3) +DS (X)*US(I,3) 

WRITE! 6 , * )  X3,Y3,Z3 
CONTINUE 


STOP 

END 


'End  of  run. 
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This  program  reads  the  geometrical  data  from  file  rec.dat  as  set 
up  by  program  rspec.f  and  generates  the  patch  coupling  coefficients 
needed  to  form  a  linear  system  of  equations*  The  numerical  integrations 
are  performed  using  Gaussian  quadrature  rules.  Two-dimensional  integra¬ 
tions  are  estimated  by  combining  two  one-dimensional  integration 
schemes,  one  embedded  inside  the  other.  The  arbitrary  interval  formula, 
Eq.  25.4.30  from  Abramowitz  and  Stegun16  was  used. 

The  primary  program  variable  definitions  are: 

A  Test  point,  loop  index.  Integer  value.  Used 

for  counting  and  ordering  rectangles. 

C  Evaluation  loop  index.  Integer  value.  Used  for 

ordering  either  conducting  or  field  boundary 
rectangles,  but  not  both  in  a  given  Do  loop. 

V1(A, I)  Absolute  (x,y,z)  position  of  vertex  1  on  rect¬ 
angle  A.  I  ■  1,2,3  for  x,y,z  components, 

respectively. 

UN(A,I),  UP(A,I)  Unit  vector  directions  n,  p,  and  s,  respec- 

US(A,I) 

tively. 

DP(A),  DS(A)  Distance  primary  and  secondary,  Dp  and  Dg ,  for 

rectangle  A. 

NP(A),  NS(A)  Number  of  primary  and  secondary  subdivisions  of 

rectangle  A. 

MODES(C)  Number  of  modes  on  field  boundary  C. 


XT(I),  XI( I) 


Test  position  and  evaluation  position  vari¬ 
ables.  I  -  1*2,3  for  x,y,z  components  of  posi- 


TSI(I), 


TS(I), 


PHI 


GPHI(I) 


H(  I) 


tlon.  Position  is  computed  in  terms  of  local 
surface  coordinates  using  Eq.  C.5. 

W(I)  Test  point  and  weight  values  for  8  point 

Gaussian  quadrature.  Used  for  numerical  approx¬ 
imation  of  self  term  coupling  coefficients.  I  - 

1*  2 i  . . . *  8. 

V(  I)  Test  point  and  weight  values  for  3  point 

Gaussian  quadrature.  Used  for  computing  all 
coupling  coefficients  except  the  self  term 
value.  I  *  1,  2,  3. 

Complex  scalar.  Numerical  value  for  $(x,x’) 
from  Section  III.  Computed  by  subroutine  GREEN 
and  depends  only  on  wavelength  and  distance 
between  x  and  x'. 

Complex  vector.  Numerical  value  for  V$(x,x') 

from  Section  III.  Computed  by  subroutine  GREEN. 

I  -  1,2,3  for  x,y,z  components,  respectively. 

♦ 

Complex  vector.  Numerical  value  for  h[i,n]  or 

4> 

h[i,n]  from  Section  III.  Computed  by  subroutine 
TEMO  for  TEmq  modes  in  a  rectangular  waveguide. 
Depends  on  local  field  boundary  position  p,s, 
the  mode  number  M,  and  on  N  ■  1,2  for  ♦ ,♦  direc¬ 
tions,  respectively.  I  -  1,2,3  for  p,s,n  direc¬ 
tion  components. 
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ZH(  I) 


HN(  I) 


ZHN( I) 


ZF(I) 


Complex  vector.  Same  as  H(l)  except  converted 


to  x,y,z  directions  by  the  transformation: 


ZH( I)  -  H(1 )UP( I)  +  H(2)US(1)  +  H(3)UN( I) ; 
I  -  1,2,3 


Complex  vector.  Numerical  value  for  3/3n  h[i,n] 


from  Section  III.  Also  computed  by  subroutine 


TEMO. 


Complex  vector.  Related  to  HN(I)  in  the  same 


way  ZH(I)  is  related  to  H(I). 


Complex  scalar.  Used  to  compute  and  store  the 


numerical  value  for 


^(c,kl]*[a,ijl 


Complex  scalar.  Used  to  compute  and  store  the 


numerical  value  for 


f  in  *  7*[a  il]J  da 

[c,kl] 


Complex  vector.  Used  to  compute  and  store  the 


numerical  value  for 


'♦la.ij]  3n  h[c,n]  "  *  V^fa,ij)^h[c,n) 

c 
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Program  DoJob 
IMPLICIT  COMPLEX  (Z) 

INTEGER  A ,  C , NP (10), NS (10) ,MODES(2) 

REAL  VI  CIO, 3) ,UN( 10,3 ),UPC 10,3) ,US( 10,3) ,DP( 10) ,DS( 10) 

REAL  XT (3) , XI ( 3 ) ,W( 8 ) ,TSI ( 8 ) ,V(3) ,TS(3) 

DIMENSION  ZFC  3 ) ,ZFI(3) ,ZFO(3) ,ZH(3) ,ZHN(3) 

COMPLEX  GPHIC3) ,PHI,H(3) ,HN(3) 

COMMON  AA,BB,AK 

DATA  TPI  /  6.283185308  / 

DATA  SMALL  /  2.SE-07  / 

Initialize  W  and  TSI  arrays  for  8  point  Gaussian  quadrature 

TSK1 )— 0.9602898565 

TSI  (2)— 0.7966664774 

TSI(3)--0. 5255324099 

TSI  (4)— 0.1834346425 

TSI C  5 ) «-TSI ( 4 ) 

TSI (61 •-TSI ( 3 ) 

TSIC7J--TSIC2) 

TSI ( 8 ) --TSI ( 1 ) 

W(l)-0. 1012285363 
WC2J-0. 2223810345 
W(3)-0. 3137066459 
W(4)-0. 3626837834 
W( 5 ) -W( 4 ) 

W( 6 ) »W( 3 ) 

W(7)-W(2) 

W( 8 ) »W( 1 ) 

Initialize  TS  array  for  3  point  Gaussian  quadrature. 

TS(l)--0. 7745966692 

TSC2J-0.0 

TS(3)--TS(1) 

V(l)-0. 5555555556 
VC2J-0. 8888888889 
VC  3 ) -VC  1 ) 

Read  rectangle  specification  data  as  generated  by  rspec.f  . 
OPENC 10, FILE- 1 rec.dat'  , STATUS- *  old  1 ) 

REWIND! 10 ) 

READ! 10,*)  NR , NRC 
DO  10  I  •  1  ,  NR 

READC 10 , * )  (V1CI,J) ,J-1,3) 

READ (10,*)  (UNC I ,J) ,J-1 ,3) 

READC 10 , * )  C  UPC  I , J ) , J- 1 , 3 ) 

READC 10 , * )  ( US (I,J),J-1,3) 

READC 10,* )  DP(I),DS(I),NP(I),NS(I) 

CONTINUE 
CLOSE ( 10 ) 

Open  file  for  vector  coupling  coefficients. 

OPENC10, FILE- ‘green. dat' .STATUS- ' unknown ‘ ) 

REWIND! 10) 
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C  Enter  run  data. 

WRITE l  *  ,  *  )  ‘Enter  ratio:  working  wavelength  /  cutoff  wavelength.' 
REAO( 5 , * )  RATIO 
AK-TPI /RATIO 
WRITE! 10  ,  * )  RATIO 

C  Enter  nuiaber  of  TE-NO  modes  at  each  field  plane. 

NFB-MR-NRC 
WRITS! 10,*)  NFS 
DO  IS  1*1, NFB 
WRITE! 6 , * )  ‘  * 

WRITE!*,*)  ‘For  field  plane  number ‘ , I ,  ‘  ' 

WRITE!*,*)  'Enter  the  number  of  TE  modes.' 

READ! S , * )  NODES! I) 

IS  CONTINUE 

WRITE! 10 ,* )  INODES! I ) ,1-1 ,NFB) 

C  Test  point  loop  over  rectangular  boundaries. 

DO  530  A-l ,NR 

WRITE! 6 ,* )  ‘Doing  work  for  rectangle  number ' , A ,'. 

HXX-DP! A ) /FLOAT! NP! A ) ) 

HYY-DS ! A ) /FLOAT ( NS ! A ) ) 

C  Test  point  loops  over  patches  on  rectangle  A. 

DO  S30  I-l ,NP! A) 

XX-HXX* ! FLOAT! I )-0 . 5 ) 

DO  530  J-l ,NS! A) 

YY-HYY* ( FLOAT ! J ) -0 . 5 ) 

C  Compute  position  of  test  point. 

DO  20  11-1,3 

XT!II)-V1!A,II )+XX*UP! A, I I )+YY*US! A, II ) 

20  CONTINUE 

C  Evaluation  loop  over  conducting  rectangles. 

DO  170  C-l.NRC 
HX-DPIC ) /FLOAT! NP(C) ) 

HY-DS I C ) /FLOAT I  NS ! C ) ) 

F»HX*HY/4 . 0 
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C  Evaluation  loops  over  patches  on  rectangle  C. 
DO  170  K-l jNP(C) 

X-HX*  (  FLOAT!  K )  -*0 . 5  ) 

DO  170  L-l ,NS!C) 

Y- HY* ! FLOAT ! L ) -0 . 5 ) 


40 


Initialize  integration  vector. 
ZG-CHPLX! 0 . 0 ,0 . 0 ) 

ZGN-CMPLX ( 0 . 0 , 0 . 0 ) 


Test  for  self  term;  or  common  patches. 

SUH-0.0 
DO  40  11-1,3 

XII II)-V1!C,II)+X*UP!C,II)+Y*US(C,II ) 

SUN-SUM+ (XT(II)-XI(II) ) **2 
CONTINUE 

IF!  ( (A.EQ.C) .AND. I !K.EQ. I) .AND. ( L . EQ . J ) ) ) . OR . ! SUM . LT . SHALL ) )  THEN 
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IS 


Do  self  torn  Integration  using  8  by  6  point  Gaussian  scheme 
DO  90  N-1,8 
AX«HX*TSI(M)/2.0+X 
ZGI-CMPLXCO. 0,0.0) 

DO  70  N-1,B 
AY-HY*TSI ( N )/2 . 0+Y 
DO  SO  II-l  ,  3 

XI ( II ) -V1CC, I I )+AX*UP(C,II )+AY*US(C,II ) 

CONTINUE 

CALL  GREEN ( XT ,XI ,PHI ) 

ZGI-ZGI+VK  N ) *PHI 
CONTINUE 

ZG-ZG+W(M)*ZGI 

CONTINUE 

ELSE 

Do  integration  using  3  by  3  point  Gaussian  scheme. 

DO  150  M-1,3 
AX-HX*TS  (  H  )  /  2 . 0-fX 
ZGI -CMPLX ( 0 . 0 , 0 . 0 ) 

ZGNI-CMPLXIO. 0,0.0) 

DO  130  N-1,3 
AY-HY*TS(N)/2 .0+Y 
DO  110  11-1,3 

XI ( II ) -V1CC, I I )+AX*UP(C,II )+AY*US(C, II ) 

CONTINUE 

CALL  GREEN ( XT , X I , PH I ) 

CALL  GGREEN ( XT , X I , GPH I ) 

ZGI -ZGI+V ( N ) *PHI 
Z-CMPLX(0. 0,0.0) 

DO  120  11-1,3 
Z-Z+UN(C,II )*GPHI( II ) 

CONTINUE 
ZGNI -ZGNI+V ( N )  *Z 
CONTINUE 

ZG-ZG+V(M)*ZGI 
ZGN-ZGN+V ( M ) *ZGNI 
CONTINUE 
ENDIF 

ZG-F*ZG 

ZGN-F*ZGN 

WRITEdO,*)  ZG ,  ZGN 
CONTINUE 

Evaluation  loop  over  field  boundary  rectangles. 

DO  350  C- ( NRC+1 ) , NR 
AA-DP(C) 

BB-DS(C) 

HX- A A/FLOAT ( NP(C ) ) 

HY- B B /FLOAT ( NS(C ) ) 

F- HX*HY/4 . 0 


C  Evaluation  loop  over  modes. 

IC-C-NRC 

DO  350  MM-1 , MODES ( IC ) 

C  Initialize  integration  sum. 

DO  190  11-1,3 

ZFC II )-CMPLX(0. 0,0.0) 

190  CONTINUE 

C  Start  summing  loops  over  patches  on  field  boundary  C. 

DO  340  I A- 1 , NP ( C ) 

X- HX* C  FLOAT ( I A ) -0 . 5  ) 

DO  340  JA- 1 , NS ( C ) 

Y • HY* ( FLOAT ( J A ) -0 . 5 ) 

C  Initialize  patch  Integration. 

DO  200  11-1,3 

ZFO( II) -CMPLX( 0 . 0 , 0 . 0 ) 

200  CONTINUE 

C  Test  for  self  patch. 

IF( IA.EQ.C) .AND. ( (I.EQ.IA) .AND. (J.EQ. JA) ) )  THEN 

C  Do  self  term  Integration  using  8  by  8  point  Gaussian  scheme 

DO  260  H-1,8 
AX-HX*TSI(M)/2.0+X 
DO  210  11-1,3 
ZFI(II)-CMPLX(0. 0,0.0) 

210  CONTINUE 

DO  240  N-1,8 
AY-HY*TSI(N)/2.0+Y 
DO  220  11-1,3 

XI(II)-V1(C,II) +AX*UP(C , I I )+AY*US l C , I I ) 

220  CONTINUE 

C  Evaluate  the  integrand. 

CALL  GREEN ( XT , X I , PH I ) 

CALL  GGREEN( XT ,XI ,GPHI ) 

CALL  TEM0 ( AX , AY , AA , BB , RATIO , 2 ,MM , H , HN ) 

Z-CMPLX(0. 0,0.0) 

DO  230  11-1,3 
Z-Z+UN(C,II )*GPHI( II ) 

ZH(II)-H(1)*UP(C,II)+H(2)*US(C,II)+H(3)*UN(C,II ) 

ZHN ( I I ) -HN ( 1 ) *UP ( C , I I ) +HN ( 2 ) *US ( C , I I ) +HN (3)*UN(C,II) 

230  CONTINUE 

DO  235  11-1,3 

ZFI( II )-ZFI( II )+W(N)MPHI*ZHN( II )-Z*ZH( II ) ) 

235  CONTINUE 

240  CONTINUE 

DO  250  11-1,3 

ZFO( II )-ZFO( II )+W(M)*ZFI ( II ) 

250  CONTINUE 

260  CONTINUE 
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Do  Integration  using  3  by  3  point  Gaussian  scheme. 

DO  320  M-1,3 
AX-HX*TS(M)/2.0+X 
DO  270  11-1,3 
ZFI(II)-CMPLX(0. 0,0.0) 

CONTINUE 

DO  300  N-1,3 
AY-HY*TS(N)/2.0+Y 
DO  280  11-1,3 

XI(II)-V1(C,II)+AX*UP(C,II )+AY*US(C,II ) 

CONTINUE 

Evaluate  the  integrand. 

CALL  GREEN ( XT , X I , PH I ) 

CALL  GGREEN ( XT , X I , GPH I ) 

CALL  TEM0 ( AX , AY , AA ,  BB  , RATIO , 2 , MM , H , HN ) 

Z-CMPLX(0. 0,0.0) 

DO  290  11-1,3 
Z-Z+UN(C,II)*GPHI(II) 

ZH(  XX )-H(l )*UP(C,XX)+H( 2 )*US(C, XX )+H( 3)*UN(C,  XX  ) 

ZHN ( II )-HN( 1 ) *UP ( C , 1 1 ) +HN ( 2 ) *US ( C , 1 1 ) +  HN ( 3 ) *UN ( C  ,  1 1  ) 

CONTINUE 

DO  295  11-1,3 

ZPI( II )-ZFI( II )+V(N)*(PHI*ZHN( 1 1 ) -Z *ZH ( 1 1  )  ) 

CONTINUE 

CONTINUE 

DO  310  11-1,3 

ZFO( 1 1 ) -ZFO( 1 1 ) +V (M)*ZFI(IX  ) 

CONTINUE 

CONTINUE 

ENDIF 

Sum  up  values  from  each  patch  integration. 

DO  330  11-1,3 

ZF( II )-ZF( II )+F*ZFO( II ) 

CONTINUE 

CONTINUE 

WRITE (10,*)  ( ZF ( II)  ,11-1,3) 

CONTINUE 

Compute  forcing  term  for  each  test  point. 

C-NRC+1 

AA-DP(C) 

BB-DS(C) 

HX- AA / FLOAT ( NP(C ) ) 

HY-BB/FLOAT ( NS ( C ) ) 

F-HX*HY/4 . 0 

Initialize  integration  sum. 

DO  370  11-1,3 

ZF( II )-CMPLX(0 . 0,0.0) 

CONTINUE 


».vv 


C  Start  summing  loops  over  patches  on  field  boundary  C. 

DO  S20  IA-1,NP(C) 

X-HX*  C  FLOAT ( I A ) -0 . 5 ) 

DO  S20  JA-1,NS(C) 

Y-HY* ( FLOAT! JA ) -0 . 5 ) 

C  Initialize  patch  integration. 

DO  380  11-1,3 

ZFO( I 1 ) -CMPLX ( 0 . 0 , 0 . 0 ) 

380  CONTINUE 

C  Test  for  self  patch. 

I F ( ( A . EQ . C ) .AND. ((I.EQ.IA) .AND. (J.EQ.JA) ) )  THEN 

Do  self  term  integration  using  8  by  8  point  Gaussian  scheme 
DO  440  M-1,8 
AX-HX*TSI!M)/2.0+X 
DO  390  11-1,3 
2FI ( I I) -CMPLX! 0 . 0 , 0 . 0 ) 

CONTINUE 

DO  420  N-1,8 
AY-HY*TSI(N)/2.0+Y 
DO  400  11-1,3 

XI! II )-Vl(C,II )+AX*UP(C,II )+AY*USCC,II ) 

400  CONTINUE 

C  Evaluate  the  integrand. 

CALL  GREEN ! XT, XI , PHI ) 

CALL  GGREEN ( XT , XI ,GPHI ) 

CALL  TEM0 I  AX , AY , AA , BB , RAT IO , 1 , 1 , H , HN ) 

Z-CMPLX10. 0,0.0) 

DO  410  11-1,3 
Z-Z+UN ( C , I I ) *GPH I ( I I ) 

ZH(II)-H(1)*UP!C,II)+H(2) *US !C,IX)+H13)*UN(C,XI) 

ZHN!  m-HN!  1)*UP!C,II  )+HN! 2 ) *US (C , II )+HN ( 3 ) *UN ( C , I I ) 

410  CONTINUE 

DO  415  11-1,3 

ZFII  II  )-ZFI(  II  )+WIN)MPHI*ZHN(  II  )-Z*ZH!  II  )  ) 

415  CONTINUE 

420  CONTINUE 

DO  430  11-1,3 

ZFO! II )-ZFO( II )+W(M)*ZFI ! II ) 

430  CONTINUE 

440  CONTINUE 


C 


390 


ELSE 

C  Do  integration  using  3  by  3  point  Gaussian  scheme. 

DO  500  M-l ,3 
AX-HX*TS ( M ) /2 . 0+X 
DO  450  11-1,3 
ZFI1II)-CMPLX!0. 0,0.0) 

450  CONTINUE 

DO  480  N-l  ,3 
AY-HY  *TS ( N ) /2 . 0+Y 
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DO  460  11-1,3 

XIdI)-Vl(C,II )+AX*UP(C,II)+AY*US(C,II) 

460  CONTINUE 

C  Evaluate  the  Integrand. 

CALL  GREEN (XT, XI , PHI ) 

CALL  GGREEN(XT,XX , GPHI ) 

CALL  TEMO ( AX , AY , AA , BB , RATIO, 1 , 1 , H , HN ) 

Z-CMPLX(0. 0,0.0) 

DO  470  IX-1,3 
Z-Z+UN ( C , XX ) *GPHI (XI) 

ZHdI)-Hd)*UP(C,II)+H(2)*US(C,II)+H(3)*UN(C,II) 
ZHN(II)-HN(!)*UP(C,II )+HN (2)*US(C,II)+HN(3)*UN(C,II) 

470  CONTINUE 

DO  475  11-1,3 

ZFKXI)-ZFI(IX)  +V (N)*(PHI*ZHN(II)-Z*ZH(II)  ) 

475  CONTINUE 

480  CONTINUE 

DO  490  11-1,3 

ZFO(II)-Z  FO ( I X ) + V ( M ) *  Z  F I ( X  X ) 

490  CONTINUE 

S00  CONTINUE 

ENDIF 

C  Sum  up  values  of  patch  Integrations. 

DO  510  11-1,3 
ZF(II)-ZF(II) +F*ZFO( I I ) 

510  CONTINUE 

520  CONTINUE 

WRITEdO,*)  CZFdl)  ,11-1,3) 

530  CONTINUE 

C  Compute  and  store  magnetic  field  variation  on  field  rectangle  C 
DO  550  C-CNRC+1 ) , NR 
AA-DP(C) 

BB-DS(C) 

HX-AA/FLOAT ( NP( C ) ) 

HY-BB/FLOAT( NS (C ) ) 

C  Loop  over  modes. 

IC-C-NRC 

DO  550  MM- 1 , MODES ( IC ) 

C  Loops  over  points  on  rectangle  C. 

DO  550  IA-1 ,NP(C ) 

X-HX* ( FLOAT ( IA ) -0 . 5 ) 

DO  550  JA- 1 , NS ( C ) 

Y-HY* C  FLOAT ( JA ) -0 . 5 ) 

CALL  TEMO (X,Y,AA,BB, RATIO, 2 ,MM,H,HN) 

DO  540  11-1,3 

ZHdI)-Hd)*UP(C,II)  +  H(2)*US(C,IX)+H(3)*UN(C,II) 

540  CONTINUE 

WRITEdO,*)  CZHdl)  ,11-1 ,3) 

550  CONTINUE 


C  Compute  and  store  incident  magnetic  field  variation. 
C-NRC+1 
AA-DP(C) 

B8-DSCC) 

HX-AA/FLOAT(NP(C)) 

HY«BB/FLOAT(NS(C) ) 

C  Loops  over  points  on  rectangle  C. 

DO  570  IA-1,NP(C) 

X-HX* ( FLOAT ( IA)-0 .5) 

DO  570  JA-1,NS(C) 

Y-HY* ( FLOAT ( J A ) -0 . 5 ) 

CALL  TEMO ( X , Y , AA , BB , RATIO , 1 , 1 , H , HN ) 

DO  560  11-1,3 

ZH(II)-H(l)*UP(C,II) +H (2)*US(C,II)+H(3)*UN(C,II) 

560  CONTINUE 

WRITEC 10  ,*  )  (ZH(II) ,11-1,3) 

570  CONTINUE 

CLOSE ( 10 ) 

STOP  'End  of  Run.' 

END 
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SUBROUTINE  GREEN ( XT , XI , PHI ) 
REAL  XT ( 3 ) , XI ( 3 ) 

COMPLEX  PHI ,  Z 
COMMON  A , B , AK 
DATA  FPI  /12. 56637062/ 

Z— CMPLXC0.0,1 .0) 

SPAN-0.0 
DO  10  1-1,3 

SPAN-SPAN+( XTC I ) -XI ( I ) )  **2 
10  CONTINUE 

SPAN-SQRT ( SPAN ) 

Z-Z*AK*SPAN 

Z-CEXP(Z) 

PHI-Z/ ( FPI*SPAN ) 

RETURN 

END 


SUBROUTINE  GGREEN ( XT , XI , GPHI ) 
REAL  XT (3) ,  X I  (  3  ) ,  D  (  3  ) 

COMPLEX  GPHI(3) ,Z,ZJ 
COMMON  A,B,AK 
DATA  FPI  /12. 56637062/ 
ZJ-CMPLX(0. 0,1.0) 

SPSQ-0.0 
DO  10  1-1,3 
D( I ) -XT( I )-XI ( I ) 

SPSQ-SPSQ+D( I ) **2 
10  CONTINUE 

SPAN-SQRT (SPSQ) 

Z-— Z J*AK*SPAN 
Z-CEXP(Z) 

Z-Z/ ( SPSQ* FPI ) 

SPAN-1 . 0/SPAN 
ZJ-CMPLXC  SPAN , AK ) 

Z-Z*Z J 
DO  20  1-1,3 
GPHI(I)-Z*D(I) 

20  CONTINUE 
RETURN 
END 
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SUBROUTINE  TEMO (X,Y,A,B,WL,N,M,H,HN) 


C  This  subroutine  computes  the  magnetic  field  components  for  TE-MO 

C  modes  in  a  rectangular  waveguide.  The  -  mode  functions  are  normalized 

C  so  that  each  mode  carries  one  watt  of  power  above  cutoff  and  one  watt 

C  reactive  (  inductive  )  below  cutoff.  The  derivative  of  the  magnetic 

C  field  components  in  the  direction  of  the  outward  normal  are  also 

C  computed. 

C 

IMPLICIT  COMPLEX  (Z) 

COMPLEX  H( 3 ) , HN ( 3 ) 

DATA  PI ,RMOE  /  3.141592654,  376.730  / 

Z J-CMPLXC 0 . 0 , 1 . 0 ) 

AK-2.0*PI/WL 
AKSQ-AK**2 
GM  *  P I *  FLOAT ( M ) / A 
GMSQ-GM**2 
WLCO-2 . 0*A/FLOAT(M) 

IF(WL.LT.WLCO)  THEN 
ZNEW-1.0 
ARG-AKSQ-GMSQ 
ZB-SQRT(ARG) 

ELSE 

ZNEW-ZJ 
ARG-GMSQ-AKSQ 
ZB— ZJ*SQRT(ARG) 

ENDIF 

ZBC-CONJG(ZB) 

ZA-4.0*ZNEW*GMSQ/(AK*ZBC*A*B*RMOE) 

AMO-REAL CCSQRT(ZA) ) 

GMX-GM*X 

H ( 3 ) -AMO*COS ( GMX ) 

HN (3)— ZJ*ZB*H(3) 

IF(N.EQ.l)  HN (33 HN ( 3  ) 

H ( 2 ) -CMPLX ( 0 . 0 , 0 . 0 ) 

HN(2)-CMPLX(0. 0,0.0) 

H  ( 1 ) -Z J*ZB*AM0*SIN ( GMX ) /GM 
HN(1)— ZJ*ZB*H(1) 

IF(N.EQ.l)  HC 1 ) — H( 1 ) 

RETURN 
END 
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This  program  reads  the  geometrical  data  from  file  rec.dat  as  set 
up  by  program  rspec.f  and  the  patch  coupling  data  from  file  green.dat  as 
generated  by  dojob.f.  The  overdetermined  system  of  equations  is  solved 
in  two  coupled  parts.  The  type  A  equations,  equations  due  to  testing  on 
the  conducting  surface,  are  partially  solved  for  the  pulse  function 
amplitudes  while  holding  the  mode  amplitudes  fixed.  That  is,  a  square 
system  of  equations  is  built  from  Eqs.  23,  24,  and  25  of  the  form  A^x^  = 
bj,  where  bj  is  a  linear  function  of  mode  amplitudes.  Next,  the  type  B 
equations,  equations  due  to  testing  on  the  field  boundaries,  Eq.  26,  are 
partially  solved  for  the  mode  amplitudes  while  holding  the  pulse  func- 


I 


tion  amplitudes  fixed.  Here  we  have  ^^2  ”  ^2 '  w^ere  is  a  linear 
function  of  the  pulse  function  amplitudes.  The  overall  iteration  takes 
the  following  form: 

1.  Hold  the  mode  amplitudes  fixed. 

2.  Make  one  ART  update  on  the  type  A  system,  A^x^  =  bj. 

3.  Hold  the  pulse  amplitudes  fixed. 

4.  Make  one  ART  update  on  the  type  B  system,  A2X2  “  ^2* 

5.  Enforce  the  conservation  of  power  constraint  on  the  mode 
amplitudes,  Eq.  28. 

6.  Repeat  until  the  solution  vector  from  the  type  B  equations 
changes  by  less  than  2  percent. 

The  ART  update  on  the  type  A  equations  is  done  using  subroutine  ARTA. 
This  subroutine  uses  a  relaxation  factor  R  that  is  between  0  and  1.  For 
the  asymmetric  iris  data  in  Section  IV,  R  ■  1,  1/2,  1/4,  1/8,  with  10 
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iterations  at  each  value  starting  with  1.  The  update  on  the  type  B 
equations  was  done  using  subroutine  ARTB  with  no  relaxation;  i.e.,  R  ■ 


1.  Conservation  of  power  was  enforced  by  subroutine  ECPMO. 
The  primary  program  variable  definitions  are: 


A,  C 

Vl( A,I) 


UN(A,I) ,  UP(A,I) 
US(A,I) 

DP( A) ,  DS( A) 

NP( A) ,  NS(  A) 

MODES(C) 

NR 

NRC 

NFB 

RATIO 


LA,  LC 


Test  point  and  evaluation  point  indices ,  respec¬ 
tively. 

Absolute  (x,y,z)  position  of  first  vertex  on 
rectangle  A.  I  -  1,2,3  for  x,y,z  components, 
respectively. 

Unit  vector  directions  n,  p,  and  s. 

Eistance  primary  and  secondary  for  rectangle  A. 
Number  of  primary  and  secondary  subdivisions  on 
rectangle  A. 

Number  of  modes  on  field  boundary  C. 

Total  number  of  rectangles  enclosing  the  volume. 
Number  of  conducting  rectangles  out  of  the 
total. 

Number  of  field  boundary  rectangles. 

Real  value.  Ratio  of  the  working  wavelength  to 
the  cutoff  wavelength  of  the  dominant  mode  on 
field  boundary  number  1. 

Integer  values.  Linear  index,  used  for  counting 
and  ordering  patches  on  rectangles  A,  C.  Com¬ 
puted  from  a  primary  index  I  and  a  secondary 
index  J  by  the  following  sequence: 


@  CONTINUE 


H(C,M,LC,I) 


HI(L1,I) 


ZG(A,LA,C,LC) 


DO  @  I-1,NP(A) 
M“NS( A)*( 1-1) 
DO  @  J-1,NS(A) 
LA-J+M 


Complex  array.  Magnetic  field  distribution  due 
to  mode  M  at  point  LC  on  field  boundary  rect¬ 


angle  C. 


1,2,3  for  x,y,z  components. 


Computed  by  dojob.f  under  variable  name  ZH(1). 
Complex  array.  Magnetic  field  distribution  due 
to  the  incident  mode  at  point  LI  on  field  bound¬ 


ary  rectangle  number  1. 


1,2,3  for  x,y,z 


components.  Computed  by  dojob.f  under  variable 
name  ZH(I). 

Complex  array.  Coupling  from  patch  LC  on  rect¬ 
angle  C  to  point  LA  on  rectangle  A.  Its  value 
is  computed  by  dojob.f  under  the  variable  name 


ZGN(A,LA,C,LC)  Complex  array.  Coupling  from  patch  LC  to  point 
LA.  Value  computed  by  dojob.f  under  variable 
name  ZGN. 

ZF(A,LA,C ,M,I)  Complex  array.  Coupling  from  mode  M  on  field 
rectangle  C  to  point  LA  on  rectangle  A.  1  - 
1,2,3  for  x,y,z  components.  Computed  by  dojob.f 
under  variable  name  ZF(I). 
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ZFI(A,LA,I) 


ZH(C,LC,I) 


ZA(C,M) 


ZHC(C,LC,I) 


ZAC(C,M) 


DA(C,LC,I) 


DB(C,LC,I) 


Complex  array.  Coupling  from  Che  incident  or 
driving  mode  in  field  rectangle  number  1  to 
point  LA  on  rectangle  A.  I  -  1,2,3  for  x,y,z 
components.  Computed  by  dojob.f  under  variable 
name  ZF(X). 

Complex  array.  Solution  vector  for  pulse  func¬ 
tion  amplitudes.  Value  for  patch  LC  on  rect¬ 
angle  C,  I  -  1,2,3  for  H  ,  H  ,  and  H'. 

P  ®  n 

Complex  array.  Solution  vector  for  mode  ampli¬ 
tudes,  mode  M  on  field  boundary  C. 

Complex  array.  Work  space,  coupling  to  equation 
LA  from  the  pulse  amplitude  coefficient  LC  on 

rectangle  C,  I  -  1,2,3  for  coupling  H  ,  H  , 

P  s 

and  IT,  respectively. 

Complex  array.  Work  space,  coupling  to  equation 
LA  from  the  amplitude  of  mode  M  on  field  bound¬ 
ary  C. 

Real  array.  Contains  denominator  terms  needed 
by  the  ART  algorithm  for  the  type  A  equations. 
Used  to  update  ZH(C,LC,I).  Computed  by  first 
access  of  subroutine  ARTA. 

Real  array.  Contains  denominator  terms  for  type 
B  equations.  Used  to  update  ZA(C,M).  Computed 
by  first  access  of  subroutine  ARTB. 
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DO  60  I-1,NP(XC) 

M-NS(IC)*(I-1) 

DO  60  J-1,NS(IC) 

LC-J+M 

READC10,*)  (  H  ( C ,  N  , LC ,11) ,11—1,3) 
CONTINUE 


XC-NRC+1 

DO  70  1*1 ,NP( IC) 

M-NS( IC)*( I-l ) 

DO  70  J-i ,NS( IC) 

LC-J+M 

READC 10,*)  (HI(LC,II ) ,11-1 ,3) 
CONTINUE 


CLOSE ( 10 ) 

WRITEt 6 , * )  'Don*  reading  data  files. 


Initialize  unknown  vectors. 
DO  80  C ■ 1 , NRC 
DO  80  1*1 ,NP(C) 

M-NS (C ) * ( I— 1 ) 

DO  80  J-l ,NS(C) 

LC-J+M 

DO  80  11-1,3 

2H ( C , LC , I I ) -CMPLX ( 0 . 0 , 0 . 0 ) 
CONTINUE 


DO  90  C-l ,NFB 
DO  90  M-l , NODES CC) 

Z A ( C , M ) -CMPLX (  0 . 0 , 0 . 0  ) 
ZOLD( C , M ) -CMPLX ( 0 .0,0.0) 
CONTINUE 


C  Initialize  Iteration  control. 

ITMAX-10 
RESMAX-S . 0E-06 
NEQ1-0 

DO  9S  A-l ,NRC 
NEQI-NEQ1+NP( A ) *NS ( A ) 

95  CONTINUE 
NEQ2-0 

DO  100  A- ( NRC+1 ) ,NR 
NEQ2-NEQ2+NP (A)* NS (A) 

100  CONTINUE 

WRITE( 6  ,*)  ‘Number  of  complex  vector  equations  - ’ ,NEQ1 , ' +' ,NEQ2 


110  IT-1 

WRITEC  6  ,  *  )  'Enter  relaxation  factor.' 
READ (5,*)  RLAX 
120  RSUM-0.0 


Start  ART  Iteration  for  type  A  equations. 

Test  loops  over  patches  on  conducting  rectangles 

aa  n  ^  a  m  _  «  a 


DO  220  A-l, NRC 
DO  220  I-l ,NP( A ) 
MM- NS ( A ) * ( I— 1 ) 
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DO  220  J-1,NS(A) 
LA* J+MM 


C  Compute  primary  component  LHS. 

ZLHS-CMPLX ( 0 . 0 , 0 . 0 ) 

DO  130  11-1,3 

ZLHS-ZLHS+UP( A,II )*ZFI (A, LA, XI ) 

DO  130  C- 1 , NFB 
DO  130  NH- 1 , MODES  CC ) 

ZLHS-ZLHS+ZA ( C , NM )*UP(A,IX ) *ZF ( A , LA ,C , NM ,11) 

130  CONTINUE 

C  Compute  coefficients  coupling  to  the  primary  term. 
DO  150  C ■ 1 , NRC 
Fl-0.0 
F2-0.0 
F3-0.0 

DO  140  11-1,3 
FI -F1+UP( A,II)*UP(C,II) 

F2-F2+UP(A,II )*US(C,IX) 

F3-F3+UPC  A , I I ) *UN (C , I I ) 

140  CONTINUE 

DO  150  K- 1 , NP( C ) 

H-NS  (C)M K-l ) 

DO  150  L- 1 , NS ( C ) 

LC-L+M 

ZHC(C,LC,1 )-Fl*ZGN(A,LA,C,LC) 

ZHC ( C , LC , 2 ) -F2*ZGN ( A , LA ,C , LC ) 

ZHC ( C , LC , 3 ) --F3*ZG  C  A , LA ,C , LC ) 

150  CONTINUE 

ZHC C A, LA, 1 ) -ZHC (A, LA, 1 ) +CMPLX ( 0 . 5 , 0 . 0 ) 

C  Do  update  to  vector  ZH. 

CALL  ARTA C A , LA , 1 , IT ,RSUM ,RLAX , ZLHS ) 

C  Compute  secondary  component  LHS. 

ZLHS-CMPLX ( 0 . 0 , 0 . 0 ) 

DO  160  11-1,3 

ZLHS-ZLHS+US (A,II)*ZFI(A, LA ,11) 

DO  160  C- 1 , NFB 
DO  160  NM- 1 , MODES ( C ) 

ZLHS- ZLHS+Z A ( C , NM )*US(A,II ) *ZF ( A , LA , C , NM ,11) 

160  CONTINUE 

C  Compute  coefficients  coupling  to  the  secondary  tei 
DO  180  C- 1 , NRC 
Fl-0.0 
F2-0.0 
F3-0.0 

DO  170  11-1,3 
F1-F1+US(A,II)*UP(C,II) 

F2-F2+US(A,II)*US(C,II) 

F3-F3+US(A,II)*UN(C,II) 

170  CONTINUE 

DO  180  K- 1 , NP ( C ) 

M-NS ( C ) * ( K-l ) 

DO  180  L- 1 , NS ( C ) 
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LC-L+M 

ZHC(C,LC, 1 )-Fl*ZGN( A ,  LA , C ,LC) 
ZHC(C,LC,2)-F2*ZGN(A,LA,C,LC) 
ZHC(C,LC,3)--F3*ZG(A,LA,C,LC) 

CONTINUE 

ZHC ( A , LA , 2 ) - ZHC ( A , LA , 2 ) +CMPLX ( 0 . 5 , 0 . 0 ) 


Do  update  to  vector  ZH. 

CALL  ARTA ( A , LA , 2 , IT , RSUM ,RLAX , ZLHS ) 


Compute  normal  component  LHS . 

ZLHS-CMPLX(0. 0,0.0) 

DO  190  11-1,3 

ZLHS-ZLHS+UN (A,II)*ZFI(A ,  LA ,11) 

DO  190  C- 1 , NFB 
DO  190  NM- 1 , NODES ( C ) 

ZLHS-ZLHS+ZA ( C , NM ) *  UN ( A , 1 1 ) *ZF ( A , LA ,C , NM ,11) 
CONTINUE 


Compute  coefficients  coupling  to  the  normal  term. 

DO  210  C - 1 , NRC 

Fl-0.0 

F2-0.0 

F3-0.0 

DO  200  11-1,3 

FI -Fl+UN ( A , I I ) *UP ( C , I I ) 

F2-F2+UN (A,II)*US(C,II ) 

F3-F3+UN (A,II)*UN(C,II) 

CONTINUE 

DO  210  K-l ,NP(C) 

M-NS ( C ) * ( K-l ) 

DO  210  L« 1 , NS ( C ) 

LC-L+M 

ZHC ( C , LC , 1 ) -F1*ZGN ( A , LA ,C , LC ) 

ZHC  C  C , LC ,2)-F2*ZGN(A,LA,C, LC ) 

ZHC  C  C , LC , 3 ) --F3*ZG ( A , LA ,C , LC ) 

CONTINUE 


Do  update  to  vector  ZH. 

CALL  ARTA ( A , LA , 3 , IT ,RSUM ,RLAX , ZLHS ) 


CONTINUE 

RA1-RSUM/FLOAT1 3*NEQ1 ) 


Start  ART  iteration  for  type  B  equations. 
Test  loops  over  patches  on  field  boundaries , 
RSUH-0.0 

DO  300  A- ( NRC+1 ) , NR 

IA-A-NRC 

DO  300  I - 1 ,NP( A ) 

MM-NS ( A ) * ( I— 1 ) 

DO  300  J- 1 , NS ( A ) 

LA-J+MH 


Compute  X  component  LHS. 

ZLHS-ZFI ( A , LA , 1 ) 

IF(IA.EQ. 1)  ZLHS-ZLHS-(0.5)*HI( LA , 1 ) 
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DO  240  C-1,NRC 
DO  240  K-1,NP(C) 

M-NS(C)*(K-1 ) 

DO  240  L-1,NS(C) 

LC-L+M 

ZLHS- ZLHS+UN ( C ,  1 ) *ZH ( C , LC ,  3 )  *  ZG( A , LA ,C, LC ) 
ZLHS-ZLHS-UP(C,1 )»ZH(C,LC,1 ) *ZGN ( A , LA ,C , LC ) 
ZLHS-ZLHS-US (C,1)*ZH(C, LC , 2 ) *  ZGN ( A , LA ,C , LC ) 
CONTINUE 

Compute  coefficients  coupling  to  the  X  component 
DO  250  C-1,NFB 
DO  250  NM- 1 , MODES (C ) 

ZAC ( C , NM ) --ZF ( A , LA  ,C , NM , 1 ) 

CONTINUE 

DO  255  NM- 1 , MODES ( I A ) 

ZAC(  I A  ,  NM  )  "ZAC  ( I A  ,  NM  )  -*-  (  0 . 5  )  *H  (  I A ,  NM  ,  LA ,  1 ) 
CONTINUE 

Do  update  to  vector  ZA. 

CALL  ARTB( I A , LA , 1 , IT ,RSUM , ZLHS ) 

Compute  Y  component  LHS. 

ZLHS-ZFI(A,LA,2) 

IF(IA.EQ.l)  Z LHS "Z LHS- (0.5)*HI( LA ,2) 

DO  260  C  - 1 , NRC 
DO  260  K-1,NP(C) 

M*NS ( C ) * ( K-l ) 

DO  260  L- 1  , NS ( C ) 

LC-L+M 

ZLHS-ZLHS+UN (CJ21*ZH(C, LC , 3 ) *ZG( A , LA ,C, LC } 
ZLHS-ZLHS-UP(C , 2 ) *ZH ( C , LC , 1 ) *ZGN ( A , LA , C , LC ) 
ZLHS-ZLHS— US ( C , 2 ) *ZH ( C , LC , 2 ) *ZGN ( A ,  LA  ,C  ,  LC ) 
CONTINUE 

Compute  coefficients  coupling  to  the  Y  component 
DO  270  C-l.NFB 
DO  270  NM- 1 , MODES ( C ) 

ZAC(C,NM)--ZF( A,LA,C,NM,2) 

CONTINUE 

DO  275  NM- 1 , MODES ( I A ) 

ZAC(IA,NM)-ZAC(IA,NM)+(0. 5 ) *H ( I A , NM , LA  , 2 ) 
CONTINUE 

Do  update  to  vector  ZA. 

CALL  ARTB ( I A , LA , 2 , IT ,  RSUM , ZLHS ) 

Compute  Z  component  LHS. 

ZLHS-ZFI ( A , LA , 3 ) 

IF(IA.EQ.l)  ZLHS-ZLHS- (0.5)*HI{LA,3) 

DO  280  C- 1 , NRC 
DO  280  K-1,NP(C) 

M-NS(C)*(K-i ) 

DO  280  L-1,NS(C) 

LC-L+M 

ZLHS-ZLHS+UN ( C, 3 ) *ZH (C,LC, 3 )*ZG( A, LA,C,LC) 
ZLHS-ZLHS-UP ( C  ,  3  )  *ZH ( C , LC , 1 ) *  ZGN ( A , LA ,  C , LC ) 
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ZLHS-ZLHS-US ( C , 3 ) *ZH ( C , LC , 2 ) *  ZGN ( A , LA , C , LC ) 

CONTINUE 

Compute  coefficients  coupling  to  the  Z  component. 

DO  290  C-l ,NFB 
DO  290  NM- 1 , MODES ( C ) 

Z AC ( C , NM ) — ZF ( A , LA , C , NM , 3 ) 

CONTINUE 

DO  295  NM-1 , MODES! IA) 

ZAC ( IA,NM)-ZAC(IA,NM)+(0.5)*H(IA,NM,LA,3) 

CONTINUE 

Do  update  to  vector  ZA. 

CALL  ARTB ( I A , LA , 3 , IT , RSUH , ZLHS ) 

CONTINUE 

RA2-RSUM/FLOAT! 3*NEQ2) 

RA3" ( RA1+RA2 ) /2 . 0 

Enforce  conservation  of  power. 

CALL  ECPM0( RATIO, DP) 

Compute  fraction  change  in  solution. 

SNUM-0.0 

SDEN-0.0 

DO  305  C-1,NFB 

DO  305  I-l , MODES (C) 

SDEN-SDEN+CABS(ZA(C,I ) ) 

ZDIFF-ZOLD ( C , I ) -ZA ( C , I ) 

SNUM-SNUM+CABS  C  ZDIFF ) 

ZOLD(C,  I  )-ZA(C, I ) 

CONTINUE 

SNUM-SNUM/SDEN 

WRITE ( 6 , * )  IT , ‘ :  Ave  error  - ' ,RA3 , ’  Soln  Chng  -*,SNUM 

IFtRAl .LE.RESMAX)  GOTO  320 
IF ( IT . GE . ITMAX )  GOTO  310 
IT-IT+1 
GOTO  120 

WRITE ( 6  ,  * )  'Hit  ITMAX  limit;  for  more  Its  type  "l",  otherwise  "O' 

READ (5,*)  NANS 

IF ( NANS .  EQ .  1 )  GOTO  110 

OPEN ( 10 , FILE- ' domore . dat ' , STATUS- ' unknown ' ) 

WRITEdO,*)  '  ' 

WRITEdO,*)  RATIO 

WRITEdO,*)  '  EA  •  '  ,RA1  ,  '  EB  -,,RA2 
DO  330  I-l  ,NFB 
DO  330  J- 1 , MODES ( I ) 

WRITEdO,*)  I,J,ZA(I,J) 

CONTINUE 
CLOSE! 10 ) 

STOP  ' End  of  Run. * 

END 
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SUBROUTINE  ARTA ( A , LA , N , ITEST ,RSUM , FAC , ZLHS ) 
IMPLICIT  COMPLEX  (Z) 

INTEGER  A ,C 

COMMON  /  B2  /  Z  H  (  8 , 4  0 , 3  ) ,  Z HC (8,40,3), DA (8,40,3) 
COMMON  /  B4  /  NP(IO) ,NS(10) , NRC , NFB , MOOES ( 2 ) 

C  Check  to  see  if  this  is  a  first  pass  access. 

IF( ITEST.GE. 2 )  GOTO  20 

C  Compute  denominator  term. 

DA(A,LA,N)-0.0 
DO  10  C-1,NRC 
DO  10  K-1,NP(C) 

M-NS(C)*(K-1 ) 

DO  10  L-1,NS(C) 

LC-L+M 

DO  10  11-1,3 
Z-ZHC ( C , LC ,11) 

DA ( A , LA , N ) -DA ( A , LA , N } +REAL ( Z  *CON JG ( Z ) ) 

10  CONTINUE 

C  Compute  the  real  part  of  the  residual. 

20  ALFA-0.0 

DO  30  C- 1 , NRC 
DO  30  K- 1 , NP ( C ) 

M-NS(C)*(K-1 ) 

DO  30  L-1,NS(C) 

LC-L+M 

DO  30  11-1,3 
ZA-ZHCCC ,LC , II ) 

ZY-ZH(C,LC , II ) 

RA-REAL ( ZA ) 

AA-AIMAG(ZA) 

RY-REAL(ZY) 

AY-AIMAG(ZY) 

ALFA- ALFA+RA*RY-AA* AY 
30  CONTINUE 

ALFA- ALFA-REAL  C  ZLHS ) 

RSUM-RSUM+ABS ( ALFA ) 

C  Update  the  vector  ZH. 

ALFA- FAC* ALFA/DA ( A , LA , N ) 

ZA-ZHC(A,LA,N) 

ZY-ZH ( A , LA , N ) 

RA-REAL (ZA) 

AA-AIMAG(ZA) 

RY-REAL ( ZY ) 

AY-AIMAG(ZY) 

RY-RY-ALFA*RA 
AY- AY+ALFA*AA 
ZH ( A , LA ,N)«CMPLX( RY , AY ) 

C  Compute  the  imaginary  part  of  the  residual. 
BETA-0.0 
DO  50  C- 1 , NRC 
DO  50  K- 1 , NP ( C ) 

M-NSCOMK-l  ) 
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DO  50  L-l jNS(C) 

LC-L+M 

DO  50  II-l ,3 
ZA-ZHC(C,LC,II) 
ZY-ZH(C,LC , II ) 
RA-REAL(ZA) 

AA-AIMAG(ZA) 

RY-REAL(ZY) 

AY-AIMAG(ZY) 

BETA- BETA+AA^RY+RA^AY 
CONTINUE 

BETA- BET A-AIHAG ( ZLHS ) 
RSUM-RSUM+ABS ( BETA ) 

Update  the  vector  ZH. 
BETA- FAC* BETA/DA C A , LA , N ) 
ZA-ZHC ( A ,  LA , N ) 

ZY-ZH ( A  ,LA , N ) 

RA-REAL(ZA) 

AA-AIMAG(ZA) 

RY-REAL(ZY) 

AY-AIMAG(ZY) 

RY«RY-BETA*AA 

AY«AY-BETA*RA 

ZH ( A , LA , N ) -CMPLX (RY , AY ) 

RETURN 

END 


SUBROUTINE  ARTBC IA,LA,N,ITEST,RSUM,ZLHS) 
IMPLICIT  COMPLEX  (Z) 

INTEGER  C 

COMMON  /  B3  /  ZA(2,8 ) ,ZAC(2 ,8 ) ,DB(2,40,3) 
COMMON  /  B4  /  NP(10) ,NS(10) , NRC , NFB , MODES ( 2 ) 
DATA  FAC  /  1.0  / 

C  Check  to  see  if  this  is  a  first  pass  access. 
IF( ITEST.GE. 2 )  GOTO  20 

C  Compute  denominator  term. 

DB(IA,LA,N)-0.0 
DO  10  C-l ,  NFB 
DO  10  NM- 1 , MODES ( C ) 

Z-ZAC(C,NH) 

DB( IA ,LA ,N ) -DB( I A , LA , N ) +REAL ( Z*CONJG( Z ) ) 

10  CONTINUE 

C  Compute  the  real  part  of  the  residual. 

20  ALFA-0.0 

DO  30  C-l, NFB 
DO  30  NM-  1 , MODES  ( C ) 

ZB-ZAC(C,NM) 

ZY-ZA(C,NM) 

RA-REAL(ZB) 

AA-AIMAGCZB) 

RY-REALCZY) 

AY-AIMAG(ZY) 

ALFA-ALFA+RA*RY-AA*AY 
30  CONTINUE 

ALFA- ALFA-REAL ( ZLHS ) 

RSUM-RSUM+ABS ( ALFA ) 

C  Update  the  vector  ZA. 

ALFA- FAC* ALFA/DB ( I A , LA , N ) 

C-  IA 

DO  40  NM-1 , MODES(C) 

ZB-ZAC(C ,NM ) 

ZY-ZA(C ,NM ) 

RA-REAL(ZB) 

AA-AIMAG(ZB) 

RY-REAL(ZY) 

AY-AIMAG(ZY) 

RY»RY-ALFA*RA 
AY-AY+ALFA*AA 
Z A ( C , NM ) -CMPLX ( RY , AY ) 

40  CONTINUE 

C  Compute  the  imaginary  part  of  the  residual. 
BETA- 0.0 
DO  SO  C- 1 , NFB 
DO  SO  NM- 1 j MODES ( C ) 

ZB-ZAC(C.NM) 

ZY-ZA(C,NM) 

RA-REAL(ZB) 

AA-AIMAG(ZB) 

RY-REAL(ZY) 


AY-AIMAG(ZY) 

BETA-BETA+AA*RY+RA*AY 

CONTINUE 

BETA-BETA-AIMAGtZLHS) 
RSUM-RSUM+ABS ( BETA ) 

Update  tha  vector  ZA. 
BETA* FAC* BETA/DB ( I A , LA , N ) 
C-IA 

DO  60  NM-l,MODES(C) 
ZB-ZAC(C,NM) 

ZY*ZA(C  jNN) 

RA-REAL(ZB) 

AA-AIMAG(ZB) 

RY-REAL(ZY) 

AY-AINAGCZY) 

RY«RY-BETA*AA 

AY«AY-BETA*RA 

ZA ( C , NM ) "CHPLX (RY ,AY) 

CONTINUE 

RETURN 

END 


I 
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SUBROUTINE  ECPMO ( WL , DP ) 

IMPLICIT  COMPLEX  (Z) 

DIMENSION  DP( 10 ) ,NM( 2 ) 

COMMON  /  B3  /  ZA(2 ,8 ) ,ZAC(2 ,8) ,DB(2,40,3) 

COMMON  /  B4  /  NP (10), NS (10), NRC ,  NFB , MODES  C  2 ) 

C  find  out  how  many  modos  can  propagate  at  aach  field  plane. 
DO  10  I-1,NFB 
IA-NRC+I 
A-DP(IA) 

NM(I)-IF1X(2.0*A/WL) 

17 ( NM ( I ) . GT . MODES ( I ) )  NM( I ) ■ MODES C I ) 

10  CONTINUE 

C  Compute  the  adjustment  factor. 

FAC-0.0 
DO  20  I-l ,  NFB 
IF(NMd).EQ.O)  GOTO  20 
DO  20  J-l ,  NM( I ) 

Z-ZA(1,J) 

FAC-FAC+REAL ( Z*CON JG ( Z ) ) 

20  CONTINUE 

FAC -SQRT (FAC) 

C  Adjust  the  solution  vector  to  force  conservation  of  power. 
DO  30  1-1 , NFB 
IF( NM( I ) . EQ. 0 )  GOTO  30 
DO  30  J-1,NM(I) 

ZA( I , J ) -ZA( I , J )/FAC 
30  CONTINUE 


RETURN 

END 
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